Perturbations of slowly rotating black holes: 
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We discuss a general method to study linear perturbations of slowly rotating black holes which 
is valid for any perturbation field, and particularly advantageous when the field equations are not 
separable. As an illustration of the method we investigate massive vector (Proca) perturbations 
in the Kerr metric, which do not appear to be separable in the standard Teukolsky formalism. 
Working in a perturbative scheme, we discuss two important effects induced by rotation: a Zeeman- 
like shift of nonaxisymmetric quasinormal modes and bound states with different azimuthal number 
m, and the coupling between axial and polar modes with different multipolar index I. We explicitly 
compute the perturbation equations up to second order in rotation, but in principle the method can 
be extended to any order. Working at first order in rotation we show that polar and axial Proca 
modes can be computed by solving two decoupled sets of equations, and we derive a single master 
equation describing axial perturbations of spin s = and s = ±f . By extending the calculation 
to second order we can study the superradiant regime of Proca perturbations in a self-consistent 
way. For the first time we show that Proca fields around Kerr black holes exhibit a superradiant 
instability, which is significantly stronger than for massive scalar fields. Because of this instability, 
astrophysical observations of spinning black holes provide the tightest upper limit on the mass of 
the photon: ^ 4 x fO"^'' eV under our most conservative assumptions. Spin measurements for 
the largest black holes could reduce this bound to m-y < 10~^^ eV or lower. 

PACS numbers: 04.40.Dg, 04.62. +v, 95.30.Sf 



I. INTRODUCTION 

Linear perturbations of black holes (BHs) play a ma- 
jor role in physics. Many astrophysical processes can be 
modeled as small deviations from analytically known BH 
backgrounds: for instance, perturbative calculations of 
quasinormal modes (QNMs) are useful to describe the 
late stages of compact binary mergers or gravitational 
collapse and BH perturbation theory provides an 

accurate description of the general relativistic dynam- 
ics of extreme mass-ratio inspirals 043 ■ Even when the 
understanding of complex astrophysical phenomena re- 
quires numerical relativity, the set up of the numerical 
simulations and the interpretation of the results arc eas- 
ier when we can rely on prior perturbative knowledge of 
the problem (see e.g. [ol. Il0|). 

Besides their interest in modeling gravitational-wave 
sources for present and future detectors in Einstein's the- 
ory and in various proposed extensions of general relativ- 
ity, perturbative studies of BH dynamics are also relevant 
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in the context of high-energy physics [llj . Perturbation 
theory can shed light on several open issues, such as the 
stability properties of BH spacetimes in higher dimen- 
sions and in asymptotically anti-de Sitter spacetimes. 
Within the gauge-gravity duality, some of the correla- 
tion functions and transport coefficients are related to 
the lowest order BH QNMs. In a semiclassical treatment 
of BH evaporation, the calculation of greybody factors 
(which may be of direct interest for ongoing experiments) 
relies heavily on our ability to understand wave scatter- 
ing in rotating BH spacetimes. 

BH perturbation theory is a useful tool to investi- 
gate issues in astrophysics and high-energy physics as 
long as the radial and angular parts of the perturbation 
equations are separable. This usually happens when the 
background spacetime has special symmetries. If sep- 
arable, the perturbation equations in the frequency do- 
main reduce to a system of ordinary differential equations 
(ODEs) • Separability is the norm if the background 
spacetime is spherically symmetric. Teukolsky [13| dis- 
covered that a large class of perturbation equations is ex- 
ceptionally separable in the Kerr metric, the underlying 
reason being that the Kerr spacetime is of type D in the 
Petrov classification. Separability is much more difficult 
to achieve in the Kerr-Newman metric [l2j and in higher 
spacetime dimensions [l^ . In the Kerr background, per- 
turbations induced by massless fields with integer and 
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half-i nteg er spins (including Dirac and Rarita-Schwinger 
fields jl2llT3l . ll5j ) are all separable, but this does not seem 
to be possible for some classes of perturbations, such as 
massive vector (Proca) perturbations p^ . 

Here we discuss a general method to study linear per- 
turbations of slowly rotating BH backgrounds that is par- 
ticularly useful when the perturbation variables are not 
separable. The method is an extension of Kojima's work 
on perturbations of slowly rotating neutron stars p7l - [l9| . 
Slowly rotating backgrounds are "close enough" to spher- 
ical symmetry that an approximate separation of the per- 
turbation equations in radial and angular parts becomes 
possible. The field equations can be Fourier transformed 
in time and expanded in spherical harmonics and they 
reduce, in general, to a coupled system of ODEs. This 
approach can be seen as a two-parameter perturbative 
expansion [23|, where the small parameters are the am- 
plitude of the perturbation and the angular velocity of 
the background. The method we present can in princi- 
ple be extended to any order in the rotation parameter; 
here we derive the perturbation equations explicitly up 
to second order. 

To first order in rotation, the final system of coupled 
ODEs can be simplified by dropping terms that couple 
perturbations with different values of the harmonic in- 
dices. This is a surprisingly good approximation: extend- 
ing previous arguments by Kojima [l9|, we show that the 
neglected coupling terms can affect the frequencies and 
damping times of the QNMs only at second or higher 
order in the rotation rate. We support the analytical 
argument by numerical results, showing that the modes 
computed with and without the coupling terms coincide 
to first order in the rotation rate. At second order in ro- 
tation, the coupling of perturbations with different har- 
monic indices cannot be neglected. However, a notion of 
"conserved quantum number" (. is preserved: perturba- 
tions with given parity and harmonic index i are coupled 
with perturbations with opposite parity and harmonic in- 
dices £±1, and with perturbations with the same parity 
and harmonic indices £ ±2 [2l| . 

The main limitation of the present approach is the 
slow-rotation approximation. Currently this is not a 
restriction in extensions of general relativity including 
quadratic curvature corrections, where rotating BH so- 
lutions are only known in the slow-rotation limit [22l - 
[26j . Furthermore, as we show in this paper, a slow- 
rotation approximation is sufficient to describe important 
effects, like the superradiant instability of massive fields 
around rotating BHs. While a second-order expansion is 
needed to describe superradiancc consistently, even the 
first-order approximation provides accurate results well 
beyond the nominal region of validity of the approxima- 
tion. Similar extrapolations have been used in the past 
to predict the existence and timescale of r-mode instabil- 
ities in the relativistic theory of stellar perturbations (see 
e.g. prl - fsoj ). and it is not unreasonable to expect that 
some quantities computed in the small-rotation regime 
(e.g. reflection coefficients) can be safely extrapolated 



to higher spin values. For all these reasons we are con- 
fident that perturbative studies of slowly rotating BHs 
will further our understanding of several interesting open 
problems in astrophysics and high-energy physics. 

As an interesting testing ground of the slow-rotation 
approximation, here we focus on massive vector (Proca) 
perturbations of slowly rotating Kerr BHs. It is well 
known that massive bosonic fields in rotating BH space- 
times can trigger superradiant instabilities |3ll437l |. For 
scalar fields the instability is well studied. It is regulated 
by the dimensionless parameter Mfi (in units G = c = 1), 
where M is the BH mass and = is the bosonic 
field masf0, and it is strongest for maximally spinning 
BHs. when Mfi ^ 1. For a solar mass BH and a field 
of mass nis ~ 1 cV the parameter Mfi ^ 10^°. In 
this case the instability is exponentially suppressed [38 1 
and in many cases of astrophysical interest the instability 
timescale would be larger than the age of the universe. 
However, strong superradiant instabilities {Mfi ^ 1) can 
occur either for light primordial BHs which may have 
been produced in the early universe [39l - l4l| or for ul- 
tralight exotic particles found in some extensions of the 
standard model, such as the "string axiverse" scenario 
|43| . In this scenario, massive scalar fields with 
10"'^^ eV < nis < 10^^® eV could play a key role in 
cosmological models. Superradiant instabilities may al- 
low us to probe the existence of such ultralight bosonic 
fields by producing gaps in the mass-spin BH Regge spec- 
trum (42l.l43j. by modifying the inspiral dynamics of com- 
pact binaries [37, 44, 45] or by inducing a "bosenova" , i.e. 
collapse of the axion cloud (see e.g. p6l - l48l |). 

Similar instabilities are expected to occur for massive 
hidden U{1) vector fields, which are also a generic feature 
of extensions of the standard model [iol - fs^ . Massive vec- 
tor perturbations of rotating BHs are expected to induce 
a superradiant instability, but an explicit demonstration 
of this effect has been lacking. While this problem has 
been widely studied for massive scalar fields [3TI - [37j . the 
case of massive vector fields is still uncharted territory; 
incursions in the topic seem to be restricted to nonro- 
tating backgrounds |16i . i53i - i55| . The main reason is that 
the Proca equation, which describes massive vector fields, 
does not seem to be separable in the Kerr background. In 
the slow-rotation approximation we can reduce the prob- 
lem of finding QNMs to a tractable system of coupled 
ODEs, where polar perturbations with angular index £ 
are generically coupled to axial perturbations with index 
£±1 (and viceversa) . 

In this paper we derive the Proca perturbation equa- 
tions up to second order in rotation, and for the first time 
we show that rotating BHs are indeed unstable to massive 
vector perturbations in the superradiant regime. At first 
order in rotation, perturbations with a given parity and 



^ In this paper, with a slight abuse of notation, we will use fj, for 
both the scalar field mass [fi = ms/K) and the vector field mass 
(/i = rriTi/h). The meaning should be clear from the context. 
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FIG. 1. (color online) Contour plots in the BH Regge 
plane [l^l corresponding to an instability timescale shorter 
than a typical accretion timescale, rsaipoter = 4.5 x 10^ yr, for 
different values of the vector field mass = ^.h (from left 



to right: 



= lO^'^'eV, 10"'-' eV, 10-^"eV, 2 x lO'^'eV). 



For polar modes we consider the S — —1 polarization, which 
provides the strongest instability, and we use both Eq. HIUOI) 
(fit II, top panel) and Eq. (|95p (fit I, middle panel), and for 
axial modes we use Eq. H95[l (bottom panel). Dashed lines 
bracket our estimated numerical errors. The experimental 
points (with error bars) refer to the mass and spin estimates 
of supermassive BHs listed in Table 2 of [S^; the rightmost 
point corresponds to the supermassive BH in Fairall 9 [53|. 
Supermassive BHs lying above each of these curves would be 
unstable on an observable timescale, and therefore they ex- 
clude a whole range of Proca field masses. 



angular index i are coupled to perturbations with i ±1 
and opposite parity (cf. Eqs. (|59l) . (|60)) and (|6T|) below). 
However we show both analytically and numerically that 
this coupling does not affect the eigenfrequencies of the 
BH. In addition, the axial sector is described by a sin- 
gle master equation [Eq. (|66p below], which is valid for 
both massive scalar and axial Proca perturbations. At 
second order in rotation the structure of the equations 
remain the same, but the £ ± 1 couplings do affect the 
eigenfrequency spectrum. 

Quasinormal frequencies, damping times of stable 
modes and instability timescales of unstable modes in 
the slow-rotation limit can be computed using standard 

methods i, [n, m, m. As a test of our approxima- 
tion scheme we revisit two problems for which "exact" 
solutions are known in the Kerr background. We study 
QNMs of massless vector perturbations to first order in 
rotation and we find that our slow-rotation method re- 
produces known results Q with accuracy better than 1% 
when the dimcnsionless Kerr parameter a = J /M"^ < 0.3. 



In addition we compute stable and unstable hound-state 
modes (which are spatially localized within the vicinity 
of the BH and decay exponentially at spatial infinity [35| ) 
of massive scalar perturbations up to second order in ro- 
tation, and we compare our results to the numerical cal- 
culations by Dolan [3^ . We find good quantitative agree- 
ment with numerical results for nonsuperradiant frequen- 
cies, and we reproduce the imaginary part of supcrradi- 
antly unstable modes within a factor 3 for a < 0.7, i.e. 
even for moderately large spin. 

Extrapolations of our perturbative results indicate 
that, because of the strong superradiant instability of 
Proca perturbations, astrophysical BH spin measure- 
ments may set the most stringent bounds on the masses 
of vector fields [s^. This is shown in Fig. [1] where we 
summarize the astrophysical implications of our results 
(cf. Sect. IVlJ for details). To be more specific, in Fig. [T] 
we show exclusion regions in the "BH Regge plane" (cf. 
Fig. 3 of (43j). i.e. we plot contours corresponding to an 
instability timescale of the order of a typical accretion 
timescale (the Salpeter time Tsaipctcr = 4.5 x 10^ yr) for 
four different masses of the Proca field (m^, = 10~^^ eV, 
10-19 eV, 10-20 eV and 2 x lO-^i eV) and for axial 
modes (bottom panel) and polar modes (top and middle 
panels). For polar modes wc consider only the S = —1 
polarization, which provides the strongest instability (cf. 
Sec. (|VI D|) for details). While our numerical results for 
the axial modes are supported by an analytical formula, 
in the polar case we have used two different functions to 
fit the numerical data at second order in the BH spin. 
The plot shows that essentially any spin measurement 
for supermassive BHs with IO'^Mq <M< IO^Mq would 
exclude a wide range of vector field masses. 

The rest of the paper is organized as follows. In Sec.HIl 
we outline the general method, which is valid for any sta- 
tionary and axisymmetric background and for any kind of 
perturbation. In Sec. lIIII we study massive Klein-Gordon 
perturbations of slowly rotating Kerr BHs up to second 
order in rotation and outline how the method can be ex- 
tended, at least in principle, to higher orders. In Sec. ITVl 
we specialize to Proca perturbations of slowly rotating 
Kerr BHs. In Sec. |V] we set up the eigenvalue problem 
for quasinormal modes and bound states, and in Sec. IVII 
wc present our numerical results. Section [VIII deals with 
the implications of the superradiant instability for astro- 
physics and fundamental physics. In Sec. IVIIII we draw 
some conclusions and we discuss possible future applica- 
tions of the method. To improve readability we relegate 
some technical material to the appendices. Appendix |X] 
lists the coefficients appearing in the scalar and Proca 
perturbation equations. Appendix [B] shows explicitly the 
perturbation equations for a Proca field on a slowly ro- 
tating Kerr background to first order in rotation, and 
appendix [Cl generalizes Dctweilcr's analytical calculation 
of the unstable massive scalar modes of a Kerr BH [g^l to 
axial perturbations of a massive vector field up to linear 
order in rotation. 

To reduce the risk of typographical errors and facilitate 
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comparison with our rcsuUs, wc made many of our calcu- 
lations available online as Mathematica notebooks [61 1. 



II. PERTURBATIONS OF A SLOWLY 
ROTATING BLACK HOLE: GENERAL METHOD 

In this section we outline the strategy to study generic 
(scalar, vector, tensor, etcetera) perturbations of any sta- 
tionary, axisymmetric background up to second order in 
the rotation parameter, but in principle the formalism 
can be extended to any order. We consider the most 
general stationary axisymmetric spacetime 

dsl = -H^dt^+Q'^dr^+r^K'^ [dd'^ + sin^ d{d(p - Ldtf] , 

where H, Q, K and L are functions of r and only. To 
second order, the metric above can be expanded as [62| 
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dsl^-F{r)[l + F2]dt' + B{t)-^ 1 + —^ dr^ 

+ r^{l + k2) [d??^ + sin^ d{dip - vodtf] , (1) 

where M is the mass of the spacetime, is a function 
of r linear in the rotation parameter, and F2, B2 and ^2 
are functions of r and "d quadratic in the rotation pa- 
rameter. Because if, K and Q all transform like scalars 
under rotation, they can be expanded in scalar spher- 
ical harmonics which, due to axisymmetry, reduce to 
the Legendre polynomials Pi{'d). As shown in Ref. j6^ . 
only t = {) and £ = 2 polynomials contribute at sec- 
ond order in rotation; therefore F2 can be expressed as 
F2(r, I?) = Fr{r) + F{){r)P2{'&), and the same applies to 
B2{r,-d) and k2{r,-d). 

At first order in rotation the background ([T]) can be 
written in the simpler form 



dsn 



B{r)-^dr^ 



F{r)di^ 
2uj{r) sin^ ddipdt 



r^d'^n 



(2) 



Note that, given a nonrotating metric, the gyromagnetic 
function vj{r) can be computed using the approach orig- 
inally developed by Hartle [62[. For example, slowly 
rotating Kerr BHs [cf. Eq. (fTS]) below] correspond to 
F{r) = B{r) = 1 - 2M/r and w = 2M'^a/r, where M 
and J = M^a are the mass and the angular momentum 
of the BH. Furthermore, the general metrics ([1]) and ([2]) 
encompass (among others) slowly rotating Kerr- Newman 
BHs and BH solutions in modified gravity theories; the 
latter have been derived analytically only up to first or- 
der [13, [IBl (and more recently to second order [26|) in 
the BH spin. 

Scalar, vector and tensor field equations in the back- 
ground metric ([T|) can be linearized in the field pertur- 
bations (that we shall schematically denote by 5X) as 
well as in the BH angular momentum. We neglect terms 
of second order in the perturbation amplitude 8X and 
terms of third order in the rotation parameter a, Fourier 



transform the perturbations, and expand them in tensor 
spherical harmonics: 



(3) 



where y^I^}!"^ is a basis of scalar, vector or tensor harmon- 
ics (depending on the tensorial nature of the perturbation 
5X) and the frequency uj is, in general, complex. The 
perturbation variables 5x'f^{r) can be classified as "po- 
lar" or "axial" depending on their behavior under parity 
transformations (i? — > tt — i?, ip — )• (yS + Tr): polar and 
axial perturbations are multiplied by (—1)^ and (— 1)''+^, 
respectively. 

With the introduction of the expansion ([3]), the linear 
response of the system is fully characterized by the quan- 
tities (5x|^(r). The perturbation equations, expanded 
in spherical harmonics and Fourier transformed in time, 
yield a coupled system of ODEs in the perturbation func- 
tions Sxfl^r). In the case of a spherically symmetric 
background, perturbations with different values of (£, to), 
as well as perturbations with opposite parity, are decou- 
pled. In a rotating, axially symmetric background, per- 
turbations with different values of m are decoupleqj but 
perturbations with different values of £ are not. However, 
in the limit of slow rotation there is a Laporte-like "se- 
lection rule" [2l| : at first order in a, perturbations with 
a given value of i are only coupled to those with ^ ± 1 and 
opposite parity, similarly to the case of rotating stars. At 
second order, perturbations with a given value of € are 
also coupled to those with ± 2 and same parity, and so 
on. 

In general, the perturbation equations can be written 
in the following form: 

Q ^ Ae, + dmAe + a^Ae 
+ a{QeVe-^i + Qe+iVe+i) 

Qe-iQeAi-2 + Qe+2Qe+iAe+2 



(4) 



0^re + dmVi + d^Vi 
+ a{QiAi^i + Qi+iAi+i) 

Qe-iQe'Pi-2 + Qi+2Qi+iPi+2 



Here 



p 



4^2 _ 1 ' 



(5) 



(6) 



Ai, Ail -^i-i -^ij -^i ^-rc linear combinations of the axial 
perturbations and of their derivatives, with multipolar 



^ From now on wc will append the relevant multipolar index £ to 
any perturbation variable but we will omit the index m, because 
in an axisymmetric background it is possible to decouple the 
perturbation equations so that all quantities have the same value 
of m. 
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index £; similarly, Vi- Ve, Ve, Vi, Vi are linear combina- 
tions of the polar perturbations and of their derivatives, 
with index i. The dimensionless parameter a keeps track 
of the order of the various terms in the slow-rotation ex- 
pansion. 

The structure of Eqs. (l4])-([5]) is the following. Pertur- 
bations with a given parity and index I are coupled to: 
(i) perturbations with opposite parity and index £ ± 1 at 
order a; (ii) perturbations with same parity and same in- 
dex (. up to order a?] (iii) perturbations with same parity 
and index ^ ± 2 at order a? . 

Prom Eq. ([6|) it follows that Q±m — 0, and therefore 
if |m| = £ the coupling of perturbations with index i to 
perturbations with indices i~ 1 and £ — 2 is suppressed. 
Since the contribution |to| = dominates the linear re- 
sponse of the system, this general property of Eqs. 
^ is reminiscent of a "propensity rule" in atomic theory, 
which states that transitions I ^ i + 1 are strongly fa- 
vored over transitions £ — )■ ^ — 1 (cf. [2l| and references 
therein) . 

Due to the coupling between different multipolar in- 
dices, the spectrum of the solutions of Eqs. (HI)-® is 
extremely rich. However, if we are interested in the char- 
acteristic modes of the slowly rotating background to first 
or to second order in a, the perturbation equations can 
be considerably simplified. We discuss the truncation to 
first and to second order in the next two sections, respec- 
tively. 



A. Eigenvalue spectrum to first order 

Let us consider the first order expansion of Eqs. (|4]) and 
We expand all quantities to first order and we ignore 
the terms Ai, Ai, Vi and Vi, which are multiplied by . 
Purthermore the terms {Vt, Ai) in Eqs. (jlj and ([5]) do 
not contribute to the eigenfrequencies at first order in 
a. This was first shown by Kojima [l9l | using symmetry 
arguments for the axisymmetric case m = 0. Here we 
extend his argument to any value of m, for the generic 
set of equations (|4]) and ([5]). Let us start by noting that, 
at first order, Eqs. (j4|) and ([5]) are invariant under the 
simultaneous transformations 

±pi-m , (7a) 
a— a, m ^ —m , (7b) 

where aim (pem) schematically denotes all of the axial 
(polar) perturbation variables with indices (£, to). The 
invariancc follows from the linearity of the terms in 
Eqs. (HI and ([5|) and from the fact that the 's are even 
functions of to. The boundary conditions that define the 
characteristic modes of the BH are also invariant under 
the transformation above (cf. Eqs. ([72|) and (|77| below). 
Therefore in the slow-rotation limit the eigenfrequencies 
can be expanded as 

u! — ojo + mujia + ujic? + 0(c?^ , (8) 



where wq is the eigenfrequency of the nonrotating space- 
time and iOn is the n-th order correction (note that Wi 
and bJ2 are generically polynomials in to but, due to the 
above symmetry, wi is an even polynomial). Crucially, 
only the terms ("P^, Ai) in Eqs. (|4]) and ([5]) can contribute 
to cji. Indeed, due to the factor a in front of all terms 
{Pii All Vi, Ae) and to their linearity, at first order in a 
we can simply take the zeroth order (in rotation) expan- 
sion of these terms. That is, to our level of approxima- 
tion the terms {Pg, Ag, Ve, Ae) in Eqs. ^ and (O only 
contain the perturbations of the nonrotating, spherically 
symmetric background. Since the latter do not explicitly 
depend on to, the m dependence in Eq. ([8]) can only arise 
from the terms {Ve, Ae) to zeroth order (recall that the 
Qf's are even functions of to, so they do not contribute). 

In the case of slowly rotating stars, this argument is re- 
inforced by numerical simulations for the to = axisym- 
metric modes [T9| and in the general nonaxisymmetric 
case [i^l , which show that the coupling terms do not af- 
fect the first order corrections to the QNMs. In this work 
we shall verify that the same result holds also for mass- 
less vector (e.g. electromagnetic) modes and for Proca 
modes of a Kerr BH. If we are interested in the modes of 
the rotating background to 0{a) the coupling terms can 
therefore be neglected, and the eigenvalue problem can 
be written in the form 

Ae + amAe = , (9) 
Ve + amVe = . (10) 

In these equations the polar and axial perturbations (as 
well as perturbations with different values of the har- 
monic indices) are decoupled from each other, and can 
be studied independently. 

B. Eigenvalue spectrum to second order 

Por what concerns the calculation of the eigenfrequen- 
cies to second order in a, the equations above can be 
simplified as follows. To begin with, we remark that all 
of the coefficients in Eqs. ^ and (O are linear com- 
binations of axial or polar perturbation functions, aem. 
Pirn, respectively. These functions can be expanded in 
the rotation parameter as 

aim - + a + a a^^ 

Plni -Pem +0- Pem +0' Pern- (H) 

The terms Ae±2 and Ve±2 are multiplied by factors a^, 
so they only depend on the zeroth-order perturbation 

functions, a^^, p^^- The terms Ai±i and Ve±i are mul- 
tiplied by factors a, so they only depend on zeroth- and 

first-order perturbation functions af2., PeH^ "^m' p'im- 

Since in the nonrotating limit axial and polar perturba- 
tions are decoupled, a possible consistent set of solutions 
of the system (H))-® has a£±2m = 0; another consis- 
tent set of solutions of the same system has p^±2m = 0- 
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Such solutions, which we can call, following Refs. [29|,|30|, 
"axial-led" and "polar-led" perturbations respectively, 
can be found by solving the following subsets of the equa- 
tions (gl)-®: 



Vi+i + amVe+i + ciQi+iAt = , 
Vii-i + amVt^i + aQiAt = , 
and 

Vi + dmVi + a^Vi + a{QtAt-i + Qt+iAi+i) = , 



(12) 
(13) 
(14) 



At+i + amAt+i + aQi+iVt = , 
Ai-i + amAi-i + aQiPi = . 



(15) 
(16) 
(17) 



In the second and third equations of the two systems 
above we have dropped the ^£±2 and ^£±2 terms, be- 
cause these only enter at zeroth order, and we have set 

"£±2m 



and p^±2>T- = 0- 



Interestingly, within this perturbative scheme a notion 
of "conserved quantum number" i is still meaningful, 
even though, for any given rotation couples terms with 
opposite parity and different multipolar index. In the 
most relevant case from the point of view of superradiant 
instabilities, i.e. the case I = m, the last equations (|14p 
and p7|) are automatically satisfied, because Qm — and 
perturbations with £ < m automatically vanish. 

It is worth stressing that, even though in principle 
there may be modes which do not belong to the classes of 
"axial-led" or "polar-led" perturbations, all solutions be- 
longing to one of these classes which fulfill the appropri- 
ate boundary conditions defining QNMs or bound states 
are also solutions of the full system (H])-® and belong 
to the eigenspectrum (up to second order in the rota- 
tion rate). This is particularly relevant since, as we shall 
show, these classes contain also superradiantly unstable 
modes. 

In the case of massive scalar perturbations, as we show 
in the next section, the coupling to £±2 can be eliminated 
by defining a suitable linear combination of the eigen- 
functions with different multipolar indices [cf. Eq. (|37p]. 
so that our procedure allows us to compute the entire 
BH spectrum up to second order. By a similar, suitable 
"rotation" in eigenfunction space it may be possible to 
recast Eqs. (H)-® in the forms (IT2|)-(Il4l) and p5|) - p7)) . 
but we did not find a general proof of this conjecture. If 
the conjecture is correct, studying the "axial-led" and 
"polar-led" systems may be sufficient to describe the en- 
tire BH spectrum at second order in rotation. 

To summarize, the eigenfrequencies (or at least a sub- 
set of the eigenfrequencies) of the general system (H))-® 
can be found, at first order in a, by solving the two de- 
coupled sets (O and ([TO)) for axial and polar perturba- 
tions, respectively. At second order in a we must solve 
cither the set (in])-([TlD or the set (IIi])-(II71) for "axial- 
led" and "polar-led" modes, respectively. We conclude 



this section by noting that this procedure can be applied 
to any slowly rotating spacetime and to any kind of per- 
turbation. For concreteness, in the next section we shall 
specialize the method to the case of massive scalar and 
vector perturbations of the Kerr metric. 



III. MASSIVE SCALAR PERTURBATIONS OF 
SLOWLY ROTATING KERR BLACK HOLES 

In order to illustrate how the slow-rotation expansion 
works, in this section we begin by working out the sim- 
plest case. We consider the massive Klein-Gordon equa- 
tion for a scalar field perturbation (jj around a rotating 
BH, and we work out the perturbation equations up to 
second order in rotation. The formalism applies to a 
generic stationary and axisymmetric background, but it 
is natural to focus on the case of interest in general rel- 
ativity, i.e. the Kerr metric in Boyer-Lindquist coordi- 
nates: 



ds 



Kerr 



= - 1- 



2Mr 



A 



d SVC? ddipdt 



(r^ +M^d^)sm^d + 



d^ sin^ ■& 



dip^ , 



(18) 



where S = -|- M'^a? cos^ A = {r — 7'+)(r — r_) and 
r± = M(l ± -\/l — a?). In what follows we shall ex- 
pand the metric and all other quantities of interest to 
second order in a. At this order, the event horizon r^^ 
the Cauchy horizon r_ and the outer ergosphere rg+ can 
be written in the form 



rs+ 



2M 1 



52 



MS 



^2 



2M 1 -cos^t?- 



(19) 



In particular, note that second-order corrections are nec- 
essary for the ergoregion to be located outside the event 
horizon. 

The massive Klein- Gordon equation reads 



(20) 



where iris = is the mass of the scalar field. We de- 
compose the field in spherical harmonics: 



E 



(21) 



and expand the square root above to second order in a. 
Schematically, we obtain the following equation: 



AiY'^ + Di cos^ = , 



(22) 



where a sum over (€, m) is implicit, and the explicit form 
of Ai and is given in Appendix I A 11 The crucial point 
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is that D( is proportional to a^, so the second term in 
the equation above is zero to first order in rotation. If 
wc consider a first-order expansion in rotation the scalar 
equation is already decoupled, and it can be cast in the 
form 



z +F— 



*f = 0, 



(23) 



where F = 1 — 2M/r, dr/dr^ = F and we defined the 
operator [l^ 



-uj^ -F 



(.{l+l) 



(24) 



Equation coincides with Teukolsky's master equa- 
tion [l3| for spin s = perturbations expanded to first or- 
der in a. The coupling to perturbations with indices £± I 
vanishes for a simple reason: Klein-Gordon perturbations 
are polar quantities, and at first order the Laporte-like 
selection rule implies that polar perturbations with in- 
dex £ should couple to axial perturbations with £ ± 1, 
but the latter are absent in the spin-0 case. At second 
order, perturbations with harmonic index £ are coupled 
to perturbations with the same parity and £±2, but this 
coupling docs not contribute to the eigcnfrcquencies for 
the reasons discussed in the previous section. Therefore 
we must solve a single scalar equation for given values of 
£ and m, that wc write schematically as 



dmVe + d^Ve ^ . 



(25) 



In order to confirm this general result, let us separate 
the angular part of Eg . (P^ . This can be achieved by 
using the identities [iTI 

cos^Y^ = Qe+iY^+^ + QfY^-^ , (26) 
sinmY' = Qe+i£Y'+^ - Qt{t + 1)Y'-^ , (27) 
cof?dY^^ {Qj+^ + Qf) Y^ 

+ Qi+iQt+2Y'+^ + QtQt-iY'-^ , 

(28) 

cosi?sin?9a^y^ = {£Q]+^ - {£ + 1)0^) Y^ 
+ Qi+iQt+2iY'^^ 

-QiQi-i{£+l)Y'-\ (29) 

as well as the orthogonality property of scalar spherical 
harmonics: 



y^y*^ dn 



The result reads, schematically, 



Ql+2Ql+lDl+2 = . 



(30) 



(31) 



By repeated use of the identity (|26p we can separate 
the perturbation equations at any order in d. Indeed, 
because of the expansion in a, only combinations of the 



form (cosi?)"y^ will appear. As we discuss in the next 
section this is also true for spin-1 or spin-2 perturbations, 
except that now the perturbation equations will contain 
combinations of vector and tensor spherical harmonics, 
and this introduces terms such as (sinz?)"9^y^, which 
can be decoupled in a similar fashion by repeated appli- 
cation of the identities listed above. This procedure is 
well known in quantum mechanics, and the coefficients 
are related to the usual Clebsch-Gordan coefficients. 
Using the explicit form of the coefficients given in Ap- 
pendix lA 11 the field equations pip schematically read 



drl 



Ue+2'^e+2 + Ue-2'^e-2 



+ Wi+2 



e+2 



dr'i 



(P'^l-2 



dri 



0, 



(32) 



where we have defined the tortoise coordinate via 
dr/dr^, = / = A/(r^ -|- a^) (expanded at second order) 
and V, U and W are some potentials, whose explicit form 
is not needed here. 

Note that the coupling to the £±2 terms is proportional 
to d^. For a calculation accurate to second order in d the 
terms in parenthesis can be evaluated at zeroth order, 
and therefore the functions ^"^^2 must be solutions of 



" ^i±2 . .(0) ^(0) 



0. 



(33) 



By substituting these relations in Eq. ([5^ we get 
d^^ 



+d^{u^%-vtlwt\) 

Finally, making use of the expressions for U and W, 
the field equations can be reduced to 



dr^ 



(35) 



where the potential is given by 
^dmulvP 



2M 



+ — — [-24A'/2 _ /^Mr {£(£ -t- 1) - 3 -f r^^) 
+2Mr^uj'^ + {£{£ + 1) + + r'^{n^ - cj^) - l) 
^r\r-2M){^Ji^-uo'-){Ql + Ql^^)] . 



(36) 



As we previously discussed, the couplings to terms with 
indices £ ±2 can be neglected in the calculation of the 
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modes. In the scalar case this can be shown exphcitly as 
foUows. If we define 



where 



2(2£- 1) 



(37) 



(38) 



then, at second order in rotation, Eq. psp can be written 
as a single equation for Zi: 



+ VtZp^ = Q, 



(39) 



which can be solved by standard methods 
coincides with Teukolsky's master equation 
s 



This equation 
13[ for spin 

perturbations expanded at second order in a. This 



is a nontrivial consistency check for the slow-rotation ex- 
pansion. In particular, the coefficients Qi in Eq. (|36p 
agree with an expansion of Teukolsky's spheroidal eigen- 
values to second order in d js^]. In fact, by extending 
our procedure to arbitrary order in d we can reconstruct 
the Teukolsky scalar potential order by order. This can 
be viewed as an independent check of the standard pro- 
cedure, which consists of expanding the angular equation 
to obtain the angular eigcnfrcqucncies (see [13] and ref- 
erences therein). 

In addition, as discussed in the previous section, the 
fact that neglecting the ^ ± 2 couplings is equivalent to a 
field redefinition implies that the entire BH spectrum of 
(massive scalar) QNMs and bound states can be found 
by neglecting those couplings. 

Finally, the near-horizon behavior of Eq. ([39|) reads 



-i knr. 



(40) 



where ku ^ - miln and Hh ^ d/{AM) + 0{d^). Note 
that, by virtue of the second order expansion, we get 
precisely ~ fc|^ close to the horizon. 

The possibility to obtain a single equation for any given 
£ and to is a special feature of scalar perturbations. The 
underlying reason is that scalar perturbations have defi- 
nite parity, so the mixing between perturbations of differ- 
ent parity cannot occur. As we show in the next section, 
this property does not necessarily hold for perturbations 
of higher spin. 



IV. MASSIVE VECTOR PERTURBATIONS OF 
SLOWLY ROTATING KERR BLACK HOLES 

The field equations of a massive vector field (also 
known as Proca's equation) read 



2 Ai^ 



0, 



(41) 



where = d^Ai,~d,jA^ and A^ is the vector potential. 
Maxwell's equations are recovered when /i = mv/h = 0, 



where is the mass of the vector field. Note that, as a 
consequence of Eq. (|1T|) . the Lorenz condition V^^*^ = 
is automatically satisfied, i.e. in the massive case there is 
no gauge freedom and the field A^ propagates 2s + 1 = 3 
degrees of freedom [l6j . 

In the ^ = case, Teukolsky showed that the equa- 
tions fo r sp in-1 perturbations around a Kerr BH arc sep- 
arable |13| . the angular part being described by spin-1 
spheroidal harmonics. In the case of a massive field the 
separation does not appear to be possible, and one is left 
with a set of coupled partial differential equations. To 
avoid these difficulties we shall consider the slow-rotation 
limit of the Kerr metric and apply the approach described 
in Sec. [ill The procedure is equivalent in spirit to that 
described for scalar perturbations, but it is more involved 
due to the spin-1 nature of the vector field. 



A. Harmonic expansion of the Proca equation on a 
slowly rotating background 

Following the notation of [H, HH we set — (t, r, x^) 
with x^ =■ We also introduce the metric of the 

two-sphere = diag(l, sin^ ■&). 

Any vector field can be decomposed in a set of vector 
spherical harmonics [6^ 

^ -(9^y^-sinl?a^yM , (42) 



sin?9 



where Y^{d,(p) are the scalar spherical harmonics. We 
expand the electromagnetic potential as follows [l6l |: 



6Af,{t,r,'d,ip) 



E 

e,m 















£.771 



"(3) 



Yi/A 



(43) 

where A = £{( + 1). Because of their transformation 
properties under parity, the functions u^^-^ belong to the 
polai^ sector when i = 1,2,3, and to the axial sector 
when i = A. In the nonrotating case the two sectors are 
decoupled The Proca equation (|^ . linearized in 

the perturbations u^^-j (i = 1,2,3,4), can be written in 
the following form: 



SUr 



A\ 



il^^ cos§ + D\ 



'B 



f cos^,?) Y' 
+ b'^P cos §] sin ddiiY'^ = , (44) 



6Iii) = (af -I- pi s\T? ■&) d^Y^ — i mPe 



sini? 



+ im + ae cos sin = , 



V — 



sini? 



= + 7£ sin^ d^Y^ + i mai 
+ (0 + Afcos?9)sinz9r^ = 0, 



2_ 

sini? 



(45) 



(46) 



where a sum over (£, to) is implicit and / denotes either 
the t component or the r component. The various radial 
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coefficients in the equations above are given in terms of radial equations: 



the perturbation functions -j in Appendix lA 21 

The perturbation equations (|44|) - (|46)) can be simpU- 
fied by using the Lorenz identity V^A^ = 0. To second Qg 
order in rotation, and for the background metric ([T5)) . 
this condition reads 



SUl = (ylf ^ + if ^ cosi? + Df^ cos^ 

+ [sf^ + i?f' cos sin ddi)Y'^ = , (47) 



where the various coefficients are again hsted in Ap- 
pendix IA21 



Each of the coefficients in Eqs. (|44p ~(|47 l) is a linear 
combination of perturbation functions with either polar 
or axial parity (cf. Appendix lA 2\ . Therefore we can 
divide them into two sets: 



D 



(I) 



£B 



(I) 



Qe+i 



D 



Qi-iQ 
+Q1+2Q1+1 



(I) 



1.-2 



• -2)B 



{I) 



Di%-{i + 3)BP, 



0, 



(50) 



Aae -imQ + Ql+^i [ipi + cj^] + Ql{£ + 1) [(£ + l)pf - a^] 
+ Qf h(^ + l)m-i - i"i((^ - l)7£-i + Af-i)] 
+ Qi+i [im+1 + im((£ + 2)je+i - Af+i)] 
-Qi-iQiii + 1) [{e ~ 2)pf_2 + 

+Qi+2Qi+it [-{t + i)pi+2 + <yi+2] = , (51) 
A/3, + i m77, + Ql^^l [£7, + \,] + Q]{i + 1) [(£ + 1)7^ - A,] 
+ Qi [-{I + + im((£ - l)p,_i + (7,_i)] 

+ Qt+i [iQ+i - im((^ + 2)pt+i - at+i)] 
-Qt-iQi{i + 1) [{t - 2)7£-2 + A,_2] 

+ Q^+2Q£+i^ [-(^ + 3)7£+2 + A,+2] = . (52) 

Note that Eqs. ([5(I l) - ([5^ have exactly the same structure 
as Eqs. (H)-®. 



Polar: A 



(i , Be , De , pe , ug, 



Axial: 77,, A« , 7,, 



where j = 0, 1, 2. In order to separate the angular vari- 
ables in Eqs. (|44|) - (|47)) we compute the following inte- 
grals: 



j6niY*^dn, {I^t,r,L); (48a) 
J 6naYt 'j-'dn , (a , 6 = ^, (p) ; (48b) 
ymaSfc*V'rff^, (a,6 = 7?,(^). (48c) 



Using the orthogonality properties of scalar and vector 
harmonics, Eq. (|30p . the relations 



Y^Y^ -/""dn = / S^Sr 7""^^^ = A(5 



(49) 



B. Proca perturbation equations at first order 

In order to make the equations more tractable, in this 
section we focus on the first-order corrections only. The 
second-order analysis is presented in Sec. lIVCl below. At 
first order, Eqs. ([50)) - ((52|) simphfy to 



4'' + Q 



A 



(I) 



+ Qe+i [A[2,-{i + 2)BP, 
Aae - i mQ 

- QfXt + i)'7£-i + Qi+i^m+i = , 

Kpt + i mrje 

- Qi{e + 1)0-1 + Qi+i£Ci+i = , 



0. 



(53) 

(54) 
(55) 



where now the coefficients are linear in a. 

To begin with, let us focus on the equations for 
monopole perturbations, £ = m = 0. The longitudi- 
nal mode of a massive vector field (unlike the massless 
case) is dynamical. Since m = 0, the monopole only ex- 
cites axisymmetric modes. In the nonrotating case these 
are described by a single equation belonging to the polar 
sector [5HI (see also [3|); for £ = 0, only the ffist two 
components u^^^ and u'^^) ^^'^ defined. However, in the 
slowly rotating case these components are coupled to the 
£ ~ 1 axial component uj^^ through Eq. ((53|) . 

When £ = m = 0, Qo = and Eq. ([53]) reduces to 



as well as the identities we find the following 



4'' + Qi a['^-2b['^ 



0, 



(56) 
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where / = 0,1,2. This is an extreme example of the 
"propensity rule" discussed in the general case: at first 
order the monopole is only coupled with axial perturba- 
tions with £ = 1. Using the explicit form of the coeffi- 
cients given in Appendix lA 21 (truncated 
at first order), the equations above can be written as a 
single equation for u^^^ : 



(f , /2(r 



drl 

2i y/^hAPujF 



MI) 



\2) 



"(4) 



(57) 



where at first order the tortoise coordinate r^, is the same 
as in the Schwarzschild case, and it is defined by dr/dr^ = 
F. The source term m^^^ is the solution of Eq. ([55]) with 

£ = 1. To first order in a we can approximate uj^j (that 
is multiplied by a in the source term) by its zeroth-order 
expansion in a, which is a solution of 



drl 



-LJ^ - F 



'(4) 



(58) 



Note that Eq. (|58|) is precisely the axial perturbation 
equation in the nonrotating case for I = \ [ig. Eqs. (|57p 
and (|58|) fully describe the dynamics of the polar t = Q 
perturbations in the slow-rotation approximation. As we 
proved in Sec. |lTl the coupling on the right-hand side of 
Eq. ([57]) does not affect the QNMs to first order in rota- 
tion. In addition, since for the monopole to = 0, to this 
order the frequency is the same as in the Schwarzschild 
case, which is extensively discussed in Ref. [T6j . 

Let us now turn to modes with ^ > 0. The equations 
for £ > at first order in a arc derived in Appendix |B] 
by using the Lorenz condition (jB2[) in order to eliminate 



The polar sector is fully described by the system 



- » 2F { 3M\ 



"(2) ~ "(3) 



r J 

A [2r^uj^ + 3F^) 



_ 2aA/Pm r 

+SF (rAFw'(2) - {r^Lu^ + AF) uf^^'j 



(59) 



V2U 



2FA 



(3) 



;:2-"(2) 



2r^uj^ul^^ + 3rF2u'(3) - 3 (A r^^) Fuf^) 



(60) 



where here and in the following a prime denotes deriva- 
tion with respect to r and the operator I?2 is defined 
in Eq. ([21)) . Note that while Eq. (pOj) only involves po- 
lar perturbations, in Eq. (|59p we also have a coupling to 

M(4) • 



On the other hand, the axial sector leads to 



'(4) 



-u 



(4) 



= -^i^M^ [{i + l)Q„„^^-i - ^Q,+i„V'+i] (61) 

(see Appendix [B] for details), where we have defined the 
polar function 



y = {A + r'^ii^) 7/(2) ~{r- 2M)u' 



(3) 



(62) 



Note that this term is similar to Eq.(25) in Ref. [l6[. As 
expected, the axial perturbation u^^-j is coupled to the 
polar functions with ^±1. Equations (|59p . (|60p and (|61l) 
describe the massive vector perturbations of a Kerr BH 
to first order in a for ^ > 0. These equations reduce to 
those in Ref. [l^ when 5, = 0. 

Since the right-hand side of Eq. (|5T|) is proportional to 
d, in our perturbative framework we can first solve for 
u^(^^ and w^^^^ to zeroth order in the rotation rate, and 
then use these solutions as a source term in Eq. (|6T|) . As 
we proved in Sec. |lTl the coupling on the right-hand side 
of Eq. (pT|) docs not affect the QNMs to first order in rota- 
tion. In Scc. lVII we shall verify this property numerically. 
Therefore, for the purpose of computing cigcnfrequcncies 
to first order in d we can simply consider the following 
polar equations: 



3A/ 
r 



\2) 



'(3) 



2dAPm 
Ar^u 



A (2r^uj^ + 3F2 



(2) 



3F [rAFu'\^) ~ {r^uj^ + AF) uf^^^ 



(63) 



2FA 



^2"(3) + 

2dA-Pm 



2r'^Lo^u\^) + irF'^u'l^) - 3 (A + r^^-j p^i_^^ 



as well as the decoupled axial equation 

AdA'Pmijj I 



(4) 



^4) 



0. 



(64) 



(65) 



Within our perturbative scheme, any cigcnfrequcncy of 
Eq. ([65]) is also a solution of the coupled Eq. (|6T|) as 
long as = ^(3)^ ~ ^' which is a trivial solution of 

Eqs. ([5^ and (pOj) for £ ± 1 at zeroth order. Further- 
more, consistently with our argument in Sec. [ill we have 
checked that there arc no other modes to order 0{d) (cf. 
Sec. EH- 

Note the similarity between Eq. (|65p . describing axial 
Proca modes, and Eq. p3|) for massive scalar pertur- 
bations. Indeed, one can show that the generalization 
of Eq. to the background metric ^ (i.e. without 
specializing to the slowly rotating Kerr metric) can be 
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written in a form that includes also massive scalar per- 
turbations. This "master equation" reads 



consistent subset of Eqs. (|50p -([52 |) reads 



FB^'I + - [B'F + F'B] 



2 2m'uj{r)uj 



\-- 


BF']^ 


)] 


\ 2r 


^ 2rFj 





= 0, 



(66) 



and it can be simplified by introducing a generalized tor- 
toise coordinate y{r) such that dr/dy = V FB. In the 
equation above s is the spin of the perturbation (s = 
for scalar perturbations and s = ±1 for vector perturba- 
tions with axial parity) . In the nonrotating case Eq. (j66p 
is exact; it also includes gravitational perturbations of a 
Schwarzschild BH if s = ±2 and F = E = 1 - 2M/r. In 
the slowly rotating case Eq. (|66|) is a complete descrip- 
tion of massive scalar perturbations for s = 0, whereas for 
s = ±1 it describes the axial sector without the £ £±1 
couplings, i.e. it is a generalization of Eq. ([55]) . 



C. Proca perturbation equations at second order 



In this section we briefly present the derivation of the 
Proca eigenvalue problem to second order in d. 

By using the Lorenz condition to eliminate the spuri- 
ous w^-^j mode, Eqs. ([5n|) - ([5^ can be written as 



2?a*a' + Va*a =0, 
I?p*P^ + Vp*p'' =0, 



(67) 
(68) 



where I?a.p are second order differential operators, Va,p 
are matrices, = (W(4)> ^(2)"^: "(3)^1 "(4)^) and *p^ = 



{u 



(2) 



/ ,/±l ,/±2 „/±2^ 
i(3)7 "(4) : "(2) > "(3) 



The explicit form of the 



equations above is quite lengthy; therefore we do not 
show it in this article, but we make it available online [6]| . 
It can be obtained using the procedure explained above 
and the coefficients listed in Appendix lA 21 The func- 
tion u^j^j can be obtained from the Lorenz condition once 
the three dynamical degrees of freedom are known. Note 
that Eqs. ([inil-dni) are particular cases of Eqs. P^ - P^ 
and Eqs. (fT5|) - pT)) . respectively. 

If we are interested in the eigenfrequencies up to sec- 
ond order in d we can drop the couplings to £ ± 2 per- 
turbations, for reasons discussed in Sec. |TT1 Therefore, a 



D 



i«, + (£-l)B«, 



+ Q,+i \Afl,-{£ + 2)B^,% 



= Aae - i mQ + Qj+i£ [£pi + en] 
+ Q^A^ + l)[i£+l)pe + ae] 
+ Qe [-{e + - im{{i - l)7£-i + A^-i)] 

+ Q1+1 [dvi+i +im{{£ + 2)ji+i - Xi+i)] , 

= Al3e + i mrje + Ql+^£ [£7^ -I- A^] 
+Ql{£ + l)[{£ + l)lt + \t\ 
+ Qi [-(£ + 1)0-1 + '^m{{£ - l)p,_i + a,_i)] 
+ Qi+i [Ki+i -im{{£ + 2)pe+i - cr^+i)] . 



(69) 



As in the scalar case, the second-order coefficients gener- 
ally contain second derivatives of the perturbation func- 
tions, i.e. ^^"(i)"^ and Since the coefficients 
are already of second order, we can use the perturbation 
equations in the nonrotating limit in order to eliminate 
these second derivatives. After some manipulation we get 
the final two sets of equations, that can be schematically 
written as 



I^kT^A' + V'aTa' = , 
P'pTp^ + V'pTp'^ = 0, 



(70) 
(71) 



where I?'a,p are second order differential operators. 



i A — 1^(4) '"(2) ' "(3) 

V'a,p are matrices. 

We remark that for a given value of £ and to, Ta and 
Tp are five- and four-dimensional vectors, respectively, 
while ^A and ^'p in Eqs. (|67p and (|68p are seven- and 
eight-dimensional vectors, respectively. This is a very 
convenient simplification, because large systems of equa- 
tions are computationally more demanding. Furthermore 
the couplings to perturbations with index £ — 1 vanish if 
£ = m, because as usual the factors = 0, and in this 
case we are left with two subsystems of dimension three. 



Tp = (U(2)'" 



''(3)''"(4) ' 



, and 



V. THE EIGENVALUE PROBLEM FOR 
QUASINORMAL MODES AND BOUND STATES 

One of the key advantages of the slow-rotation approx- 
imation is that the perturbation equations can be solved 
using well-known numerical approaches. The integration 
proceeds exactly in the same way as in the nonrotating 
case. For example, for Proca perturbations of a slowly ro- 
tating Kerr BH we can compute the characteristic modes 
by following a procedure similar to the Schwarzschild case 
discussed in Ref. [l^, as follows. 

The perturbation equations (|70p - ([7T|) . together with 
appropriate boundary conditions at the horizon and at 
infinity, form an eigenvalue problem for the frequency 
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spectrum. In the near-horizon hmit, X'a,? cP/dr^ 
and Va,p — > (o; — mfi//)^, so that at the horizon we 
impose purely ingoing-wave boundary conditions: 



where w^.j^ is a constant and 



kn ~ uj — mflff = uj h 0(a'^) . 

4M ^ ' 



(72) 



(73) 



We have introduced the horizon angular frequency iln = 
a/(2A'/r+) and expanded it to second order. 

In the massless case, due to the boundary condi- 
tion ([7^ . superradiant scattering for scalar, electromag- 
netic and gravitational perturbations is possible when 
ujji < mQ,H (66| . i.e. (to second order in rotation) when 



a > 



AMwr 



m 



(74) 



where ujr is the real part of the mode frequency, i.e. 
u! = LOfi + iu>i. Superradiance is possible because the 
energy flux at the horizon 



Er^ = lim 



(75) 



is negative when kn < 0. In the equation above T^^^ is 
the stress-energy tensor of the perturbation. The same 
argument can be applied to massive scalar perturbations. 

The case of massive vector perturbations is more in- 
volved. For purely axial perturbations close to the hori- 
zon, by using Eq. (|75p and the stress-energy tensor of a 



Proca field we get E^^ = 



Et with 



Ei 



Lok 



H 



M{e + i)M 



6l"(4)ffl 



0{d 



(76) 



which shows that the energy flux across the horizon is 
negative when kn < 0. The general case involves terms 
proportional to each component w^^^^, as well as terms 

proportional to 

As we discuss in Sec. IVII below, our results turn out 
to be very accurate for moderate values of d and they 
are reliable even in the superradiant regime, defined by 
Eq. (fM]) . as long as ujuM ^ 1. Superradiant scat- 
tering leads to instabilities for massive scalar perturba- 
tions [1^ [gOI- In the Proca case, the numerical results 
discussed in Sec. |Vl] below show that, when the super- 
radiant condition (j74p is met, the imaginary part of the 
modes crosses zero. Thus our numerical data show hard 
evidence, for the first time, that massive vector fields 
trigger (as expected) a superradiant instability. 

The asymptotic behavior of the solution at infinity 
reads 



where k^o = — (^^, so that Re[koo] > 0. The bound- 
ary conditions = yield purely outgoing waves at 
infinity, i.e. QNMs If instead we impose C(j) = we 
get states that are spatially localized within the vicinity 
of the BH and decay exponentially at spatial infinity, i.e. 
bound states (see e.g. [3, |35|). Stable QNMs are more 
challenging to compute than bound states. In the former 
case the boundary conditions above imply that purely 
outgoing waves blow up as oo, whereas purely in- 

going waves are exponentially suppressed at infinity. For 
bound states, direct integration of the equations com- 
bined with a shooting method is sufficient, but QNMs 
are more efficiently computed by other means, for exam- 
ple via continued fraction methods [3]. 



A. On the superradiant regime and the second 
order expansion 

Here we comment on some important features that 
would be missed in a first order treatment. First, at 
second order the structure of the background metric is re- 
markably different, because all metric coefficients acquire 
0{a^) corrections. This affects the location of the event 
horizon, of the ergosphere and of the inner Cauchy hori- 
zon [cf. Eq. (|19l) ]. Second-order corrections are known to 
approximate the Kerr solution much better than a first- 
order expansion (see e.g. [67|). 

From a dynamical point of view, axisymmetric modes 
acquire second-order corrections that break the m = 
degeneracy of the first-order case. Most importantly, at 
second order the superradiant regime of vector and scalar 
fields can be described within our perturbative approach 
in a self- consistent way. The superradiant condition (|74p 
may look like a first-order effect. However, at the onset 
of superradiance uiM ^ d (at least in the most relevant 
cases, i.e. when m is of order unity). In this case, terms 
like o;^ in the field equations (which are crucial, in partic- 
ular in the study of the mode spectrum) are of the same 
order of magnitude as second order quantities. 

The field equations are of second differential order, so 
the linearized field equations (at first order in a) contain 
at most terms of order 



aojM . 



(78) 



In principle, terms of order a(a;il/)^ would also be al- 
lowed, but those do not appear in the linearized Proca 
equations (|50p - (|52p . In the massless case this is consis- 
tent with a first-order expansion of the Teukolsky equa- 
tion, which does not contain any a(a;M)^ term for scalar, 
vector and tensor perturbations. 

Thus, if LoM ~ a the second and fourth terms listed 
above would be as large as second-order terms in a, which 
are neglected in a first-order expansion. In a second-order 
J expansion, instead, one also keeps terms proportional to 

, a?. This would be enough to consistently describe the 

(77) superradiant regime up to first order in a, i.e. up to terms 
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^ auj and ^ uP'. For a related discussion in the case of 
neutron star r- modes, we refer the reader to Ref. [68ll. 



B. Continued fraction method for quasinormal 
modes and bound states 

From a conceptual point of view, numerical calcula- 
tions in the slow-rotation approximation are not more 
complicated than in the nonrotating case. For example, 
in order to apply the continued fraction method we can 
write down a recurrence relation similar to Eqs. (31) and 
(32) of Ref. [H], starting from Eqs. (fSQ]) . ([60l) and (|65|) . 
by imposing the ansatz 



where U„ = (a„ , ) is a two-dimensional vectorial 
coefficient and «„, /3„, 7„, (5„, p„ and cr„ are 2x2 
matrices, whose explicit form we do not present here for 



By using a matrix- 
701 the system above 



(r - r+)-^^^"r^e-''^ ^ " , (79) 



where v = —q + lo^ Iq + 2ikH and q = for bound 

states and for QNMs, respectively. For the axial equation 
this ansatz leads to a three-term recurrence relation of 
the form 



0, 



brevity but it is available onlin e 161 
valued Gaussian elimination [6 
can be reduced to a three-term matrix-valued recurrence 
relation, which can be solved with the method discussed 
in Ref. [H]. 

The continued fraction method works very well for 
both QNMs and bound-state modes. We computed 
QNMs and checked that they yield the correct limit in 
the massless case (see Sec. IVIBI below), but in the fol- 
lowing we will focus mainly on bound states, that can 
become unstable in the superradiant regime. 



Direct integration and Breit-Wigner resonance 
method for bound states 



To compute bound state frequencies it is possible to use 
either a direct integration method or the Breit-Wigner 
resonance method. We start with a series expansion of 
the solution close to the horizon: 



whose coefficients (to first order in a, and setting M 
in these equations only) read 

-4{l + n)q^{l + n 
4q {q{l+£{£ + l)^ 
+n{2 + 6q) - s^) ~ 
-{l + 2n + l2q)LO^ 
+4i dmq (q + 2nq -\ 



-r 



(83) 



where is expanded to second order and the coefficients 



- 4iw) 

2^2 + q{3 + 4q) 
4iq{l + 2n + 3q)uj 
+ 4iw3) 

3q^ - 2i(?w - io"^) , 



4iam(l -|- n)q'' (80) 6n'' (n > 0) can be computed in terms of b^' by solving 



(81) 



7n 



-4 [{q{n + q) 
'4idmq (q{n 



,2\2 



- 2i quj 
q) — 2\quj — I 



2 21 

s q \ 



(82) 



for s = ±1. If we set s = 0, the equations above are 
also valid for massive scalar perturbations of a Kerr BH 
in the slowly rotating limit because, as we have shown 
above, there exists a single master equation describing 
both massive scalar and axial vector modes [cf. Eq. (|66p ]. 
We have investigated the massive scalar case to test the 
robustness of our results, because in this case the per- 
turbation equations on a generic Kerr background are 
separable |60j and the eigenvalues can be computed by a 
direct solution of the Teukolsky equation 0, HH . 

In the slowly rotating case the polar sector leads to a 
six-term, matrix- valued [l6j recurrence relation 

aoUi + /3oUo = , 

aiUa +/3iUi +7iUo = 0, n>0, 

aaUs + ;32U2 + 72U1 + J2U0 = , n > 1 , 

a3U4 + /33U3 + 73 U2 + h^i + P3U0 = , 71 > 2 , 
anU„+i + /3„U„ + 7„U„_i + (5„U„_2 

-l-p„U„_3 + cr„U„_4 = , n > 3 , 



the near-horizon equations order by order. In the direct 
integration method, the field equations are integrated 
outwards up to infinity, where the condition C(j) = 
in Eq. ([77]) is imposed (see Ref. [ll] for details). 

The Breit-Wigner resonance method, also known as 
the standing- wave approach [gsI . [7ll - l73j . is well suited 
to computing QNMs and bound states of the system of 
equations (|59p - (|5T|) in the case of slowly damped modes, 
i.e. those with cj/ ^ lor. In this case the eigenvalue 
problem can be solved by lookin g fo r minima of a real- 
valued function of a real variable [72| . We briefly explain 
the procedure below, where we extend it to deal with a 
system of coupled equations. For clarity we only consider 
first-order corrections, but our argument applies also to 
the second-order case described by Eqs. ([7(H) - (|7ip and, 
in principle, to any order. 

Since ^ = 0, 1, 2, .., the full system ^ and (EH) 
formally contains an infinite number of equations. In 
practice, we can truncate it at some given value of ^, 
compute the modes as explained below, and finally check 
convergence by increasing the truncation order. Let us 
suppose we truncate the axial sector at £ = L and the 
polar sector &t I ~ L + \ . i.e. for a given m we assume 



0. 



u) 







when 



£> L. 



(84) 



with J = 1,2,3 denoting the polar perturbations. All 
perturbations vanish identically for £ < |m|. 

When m = 0, the truncation above reduces the sys- 
tem to iV = 3L coupled second-order ODEs for L — 1 
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axial functions and 2L + 1 polar functions, including the 
monopole, described by Eq. ([57)) . When |m| > the 
truncated system contains N = iL — 'd\m\ + 2 second- 
order ODEs (for L— \m\ axial functions and 2L — 2\m\ + 2 
polar functions). In all cases we are left with a system 
of N second-order ODEs for N perturbation functions, 
which we collectively denote by y(p) (p = 1, ...,N). 

At the horizon each function is described by ingoing 
and outgoing waves. We impose a purely ingoing wave 
boundary condition analogous to Eq. 



Vip) 



(85) 



where again the coefficients ci^'' {n > 0) can be com- 
puted in terms of Cq^' . A family of solutions at infinity 
is then characterized by N parameters, corresponding to 
the TV-dimcnsional vector of the near-horizon coefhcients, 
Co = {c^^}. At infinity we look for exponentially decay- 
ing solutions, which correspond to bound, slowly damped 
modes. The spectrum of these modes can be obtained as 

follows. We first choose a suitable orthogonal basis for 

(v) 

the A^-dimensional space of the initial coefficients Cq 
(see also Ref. 0]). We perform N integrations from the 
horizon to infinity and construct the N x N matrix 



Sm(a;) = lim 



(2) 

2/(2) 



^(1) \ 



(86) 



where the superscripts denote a particular vector of the 



basis, i.e. 



V 



corresponds to Cq = {1, 0, 0, 0}, yfj^ 
corresponds to Cq = {0, 1,0, ...,0} and y^^^-^ corresponds 

to Co = {0, 0, 0, 1}. Finally, the bound-state frequency 
uji^) = ujR -\- iLOj IS obtained by imposing 



det S,„(i:jo) = . 



(87) 



So far we have not imposed the Brcit-Wigner assumption 
LJi <^ ujfj. In fact, the procedure discussed above can be 
used to perform a general direct integration in the case 
of coupled systems. 

By expanding Eq. (|87)) about uir and assuming loj <C 
ujR we get 



detSm(a;o) ~ det Sm(wfl) -fiw/ 



rf[dctS„(w)] 



du 



0, 



which gives a relation between d[detSm{'-^)]/duj and 
det S„i at w = cj^. We consider the function detSm 
restricted to real values of w. A Taylor expansion for 
(real) u! close to ujji yields (using the relation above): 



det S,„(w) ~ detSm(a;/?) 



1 



UJ — LJR- 



lo;/ . 
(89) 



Therefore, on the real-o; axis, close to the real part of 
the mode. 



I det S,„ (aj)p oc (oj — ujb) + lo 



I ■ 



(90) 



To summarize, to find the slowly damped modes it is suf- 
ficient to integrate the truncated system N times for real 
values of the frequency w, construct the matrix S„i (cj) 
and find the minima of the function | det Smp, which rep- 
resent the real part of the modes. Then the imaginary 
part (in modulus) of the mode can be extracted through 
a quadratic fit, as in Eq. (|90p . We note that the same 
method can be straightforwardly extended to compute 
the slowly damped modes of the general systems ([E 
(HI and dHll-dlTl). 
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No coupling 


i = 1 (DI) 


0.099484532 


1.1 • lO"'' 


Full system, L = 


a,2,3,4(BW) 


0.099484563 


1.1 • 10"' 


No coupling 


i = 1 (DI) 


0.09987320753 


4 ■ 10"" 


Full system, L 


= 2, 3, 4 (BW) 


0.09987183250 
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TABLE I. Examples of polar (top rows) and axial (bottom 
rows) bound-state modes for m — 1, Mfj, = 0.1 and a = 0.1 
computed at first order. Modes were computed via direct 
integration (DI) of the system (|59|l - (|60|l without couplings 
I ^ I ± 1, as well as with the Breit-Wigner (BW) method 
applied to the full system (|59|l - (|60|) and (|61|l . for different 
truncation orders L. The modes are insensitive to the trun- 
cation order within the quoted numerical accuracy. 

By applying the Breit-Wigner method we have numer- 
ically verified the argument given in Sec. |lTl i.e. that the 
£ — £ ± 1 couplings do not affect the eigenfrequencies 
in the slow-rotation limit. As an example, in Table |T] 
we compare two modes (m = 1, a = 0.1, M^, = 0.1) 
of the full system computed with the Breit-Wigner res- 
onance method for several truncation orders L and the 
same modes computed with a direct integration of the 
system of equations without the £ — > ^ ± 1 couplings. Up 
to numerical accuracy the mode is insensitive to the trun- 
cation order and, most importantly, it agrees very well 
with that computed for the system without couplings. 
Note that L = 1 corresponds to the uncoupled polar sys- 
tem with £ = 1. Therefore the small discrepancy is not 
due to the coupling terms, but to some inherent numer- 
ical error of the resonance method, which becomes less 
accurate when the imaginary part of the mode is tiny. 



VI. RESULTS 

A. Numerical procedure 

In principle, by integrating the full systems ([70]) and 
(|7ip for a given truncation order L and a given value of m 
one can obtain the full spectrum of quasinormal modes 
and bound modes for both the axial and polar sectors 
and for any £ < L. We have computed bound states 
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and QNMs via the Breit-Wigner procedure of Sec. IV CI 
We double- checked the results using two additional, inde- 
pendent techniques, also described above: the continued 
fraction method and direct integration of the reduced 
equations ((70)) and ([TT]) . The results agree within nu- 
merical accuracy. The Breit-Wigner procedure and the 
continued fraction method have been implemented up to 
first order in the rotation parameter; the direct integra- 
tion, instead, has been extended up to second order in 
a, in order to validate the results of the first-order inte- 
gration, as we shall discuss below. Furthermore, when 
LOi <^ ujR we found very good agreement with the Breit- 
Wigner method (cf. Table U). 

Exploring the full parameter space of the eigenfrequen- 
cies at second order in d is numerically demanding. For 
this reason, in Sec. IVI Dl below we will present a more 
extensive survey of results at first order, and compare 
them to the second-order calculation in Sec. IVI El in se- 
lected cases. 



B. A consistency check: massless vector 
perturbations 

As a preliminary test of our method, we have com- 
puted the QNMs of massless vector (i.e., electromagnetic) 
perturbations of a Kerr BH to first order in a. These re- 
sults can be compared with those obtained by solving the 
Tcukolsky equation without imposing the slow-rotation 
approximation (see e.g. Q)- In the massless limit the 
axial equation (|65|) reduces to 



10' 



(4) 



AdM^muj 



1) 



F^u(4) =0. (91) 



For the polar sector in the massless limit, we can define 
exactly the same master variable as in the nonrotating 
The polar sector in the massless limit is de- 
a fourth-order equation. As in the nonrotat- 
one solution of this equation is a pure-gauge 
In the slow-rotation ap- 



I4|- 



case 

scribed b:v 
ing case 

mode and can be eliminated, 
proximation we can recast the other solution in terms of 
a master function, which also satisfies Eq. (|?T|) . There- 
fore axial and polar perturbations have the same spectra. 
This is consistent with the fact that electromagnetic per- 
turbations of the full Kerr geometry are described by a 
single master equation in the Teukolsky formalism 

Due to isospectrality, in the massless limit we only need 
to solve Eq. (|9ip . The modes can be computed via the 
continued fraction method introduced above, where the 
coefficients of the recurrence relation can be obtained by 
setting ^ = and s = ±1 in Eqs. (|M)) - (|5^ . 

We compared our results with the exact massless vec- 
tor modes of a Kerr BH, computed by solving the Teukol- 
sky equation, for (. = 1 and different values of m (see 
Fig. [2]). The first order approximation performs very 
well, even for relatively large values of the BH rotation 
rate. Remarkably, the real and imaginary parts of the 



^ -2 



— "r' 


m=l 


-- Mj, 


m=l 


— "r' 


m=0 


-- CO,, 


m=0 


— 


m=-l 


-- CO,, 


m=-l 



1_%_error__ 




FIG. 2. (color online) Percentage difference between QNMs 
of massless vector perturbations in the slow-rotation limit at 
first order and the "full" numerical solution of the Teukolsky 
equation in the Kerr metric Q. The deviation scales like a"^ 
as a ^ 0. 



QNMs computed with our approach deviate by less than 
1% from their exact values if a < 0.3. 
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FIG. 3. (color online) Comparison between the exact, 
Teukolsky-based result and the results obtained by our slowly 
rotating approximation at first and second order for the imag- 
inary part of the scalar fundamental bound-state mode with 
I = m = 1 and ^iM = 0.1. For comparison, we have also com- 
puted the same mode as obtained by expanding the Teukolsky 
equation at third and fourth order. 



C. A second test: bound state modes for scalar 
perturbations 

We can also investigate the accuracy of the slow- 
rotation approximation for massive scalar perturbations, 
another case in which the Tcukolsky equation can be 
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solved exactly (see e.g. [35|). In Fig.|3]we show the imag- 
inary part of the bound state modes with £ = m = 1 
and Mfi ~ 0.1, computed by solving Eq. (p9|) at first 
and second order via direct integration, against the ex- 
act Teukolsky-based result obtained with the continued 
fraction method [3, [H]. For comparison, we also show 
the results obtained by solving Teukolsky's equation at 
third and fourth order. As shown in Fig.|31 the imaginary 
part crosses the axis when the superradiant condition is 
satisfied. In the stable branch (a < 0.4) even first-order 
results are in good quantitative agreement with the "ex- 
act" calculation. In the unstable branch (a > 0.4) the 
first-order approximation is in only qualitative agreement 
with the exact result, but the second-order approxima- 
tion is in quantitative agreement with the numerics at 
the onset of the instability. It is still quite remarkable 
that even first-order results can correctly reproduce the 
onset of the instability. We shall return to this consider- 
ation below, when we will deal with the massive vector 
case. 
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FIG. 4. (color online) Bound state modes obtained with a 
first-order Breit-Wigner method applied to the full system. 
We show the determinant | det Sm | as a function of the real 
part of the frequency for M/j. — 0.1 and 2 — 0.1. According 
to Eq. (|93|l , the real part of modes with the same £ + n + S is 
approximately degenerate for Mjj, <^ 1. 



^ = 0, S* = 1, in a gree ment with the rules for addition of 
angular momenta [l6| . In the limit a = we recover this 
scaling. We have also analyzed the a-dependence of the 
modes for small d. When M/x ^ 1 we find the expected 
hydrogen-like behavior, Wfl, ~ /x, which remains valid also 
in the slowly rotating case. When Mfj, < 0.1 the real part 
of the modes is roughly independent of a, and it is very 
well approximated by the relation 



2 2 



1 - 



e+n+s+i 



(93) 



where n > is the overtone number (cf. |35|). For the 
axial case (5* = 0) this relation is validated by the ana- 
lytical results presented in Appendix [C] [sec in particular 
Eq. (|C9p ]. where we solve the axial equations in the limit 
A//i <C 1. Equation (|93p is also supported by the nonro- 
tating result given in Eq. (49) of [ig (see also H^). 

Equation (j93p predicts a degeneracy for modes with 
the same value oi £ + n + S when Mn <C 1. In the Breit- 
Wigncr method, the mode frequencies can be identified 
as minima of the real- valued function | dot SmP- The de- 
generacy predicted by Eq. ((93l) is not exact for Mfi = 0.1, 
as illustrated in the inset of Fig. 21 where we display the 
minima of |detS„ip when d = 0.1. The first minimum 
on the right corresponds to£ + n + S~0, which can 
only be achieved for the fundamental polar mode with 
{i,n,S) = (1,0,-1). When^-l-n+S' = 1 we have a three- 
fold degeneracy, corresponding to {i,n,S) = (1,0,0), 
(2, 0, —1), (1, 1, —1). The approximate nature of the de- 
generacy is shown in the inset, where three distinct (al- 
beit very close) minima appear. For £ + n + S = 2 there 
is a five-fold degeneracy, which can be resolved with high 
enough resolution. Note that according to Eq. (^5]) the 
imaginary part of the modes is tiny when £+n+S is large. 
This makes it difficult to numerically resolve higher over- 
tones and modes with large £. 

The imaginary part of nonaxisymmetric modes shows 
the typical Zeeman-like splitting for different values of 
m when a ^ 0, as shown in Fig. Ofor the axial modes 
and for the polar mode with S = —1. For m > the 
imaginary part of the frequency decreases (in modulus) 
as 5, increases, and it has a zero crossing when 



D. Proca modes to first order 

Unlike the massless case, axial and polar modes for 
massive vector perturbations are not isospectral. Fur- 
thermore there exist two classes of polar modes, which 
can be distinguished by their "polarization" [16|. For 
small Mfi, Rosa and Dolan found that the imaginary 
part of the bound states of a Schwarzschild BH scales as 



uji ~ fi{Mfi) 



M+5+2S 



(92) 



where 5 = for axial modes and 5 = ±1 for the two 
classes of polar modes. The monopole corresponds to 



(94) 



which according to Eq. ([74)) corresponds to the onset of 
the superradiant regime. Recall that a second-order cal- 
culation is needed to describe the superradiant regime in 
a self-consistent way, because the latter is well beyond 
the nominal regime of validity of the first-order approx- 
imation. It is quite remarkable that even the first-order 
approximation predicts that the instability should "turn 
on" at the right point, as shown in Fig. [5] for £ ^ m ~ 1 
and several values of /i. The quantitative accuracy of the 
first-order approximation is questionable in the superra- 
diant regime. However, we will now show that second- 
order results indicate that we are correctly capturing the 
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FIG. 5. (color online) Absolute value of the imaginary part of the axial (left panel) and of the polar S = — 1 (right panel) 
vector modes as a function of the BH rotation rate a for £ — 1, Mfi — 0.05 and different values of m, computed at first order. 
For m = 1, the modes cross the axis and become unstable when the superradiance condition (|94|l is met. In the left panel, the 
red dot-dashed line denotes the analytical result (|C12|) . 



order of magnitude of the instability even for moderately 
large rotation. 

E. Proca modes to second order 

As discussed in Sec. IIV C[ the perturbation equations 
for massive vector fields at second order in rotation, 
Eqs. ([70)) and ([7T|) . are obtained from Eqs. (|B0)) . using 
the perturbation equations in the nonrotating limit to 
eliminate the second derivatives. The explicit form of 
the equations is available online [6l| . 

Solving the eigenvalue problem for these equations is 
numerically more demanding than in the case of pertur- 
bations at first order in a, especially when the imaginary 
part of the modes is tiny (as in the axial case and in the 
small- mass limit). For this reason we present only a se- 
lection of results that allow us to validate the accuracy of 
the first-order approximation, and (more in general) to 
estimate the errors due the slow-rotation approximation. 

The imaginary part of the fundamental S' = axial and 
S = — I polar modes as a function of d at first and second 
order are shown in Fig. [T] The slow-rotation method cor- 
rectly predicts the onset of the Proca instability already 
at first order, where (a priori) large deviations could be 
expected. Furthermore, the value of a which corresponds 
to the onset of the instability is slightly smaller at second 
order, consistently with the fact that the horizon location 
gets negative second-order contributions [cf. Eq. (jl9p ]. 
Fig. [7] actually displays a remarkably good quantitative 
agreement between the first- and second-order calcula- 
tions. The main difference is that in the superradiant 
region the instability predicted by the first-order expan- 
sion is weaker (by a factor of a few) than the instability 
computed at second order. Another important point is 
that our perturbative approach can consistently describe 



the superradiant regime as long as Mtu <ti 1 and d <^ 1. 
Due to the hydrogen-like spectrum of the bound states, 
this means that we are limited to considering small values 
of fiM. Again, this expectation is validated by Fig. [71 

At second order, all curves in Fig. [T) show a maximum 
at large values of the spin. This maximum is an arti- 
fact of the second-order approximation: a similar fea- 
ture appears also in the scalar case in which the exact, 
Teukolsky-based result is monotonic (cf. Fig. [3]). In what 
follows we have discarded numerical data beyond such 
maxima, and we have taken this truncation into account 
when estimating numerical errors. 

The bottom line of this discussion is that the first or- 
der calculation provides a reliable estimate of the order 
of magnitude of the instability. The second-order calcu- 
lation should be quantitatively reliable even in the unsta- 
ble regime, at least up to moderately large values of d. 
Roughly speaking, the reason why the first-order calcu- 
lation captures the correct qualitative properties of the 
instability can be traced back to the fact that the super- 
radiant instability condition uj < mVLH is a "first order 
effect" . Because is an odd function of a, a first-order 
approximation of this relation can be expected to be valid 
modulo third-order corrections. 



F. The minimum instability growth scale 

The agreement between first- and second-order calcu- 
lations can be used to estimate the order of magnitude 
of the instability by extrapolation. We will fit our data 
at the onset of the instability, and then extrapolate to 
the superradiant region. This procedure should give us 
at least a correct order of magnitude for the instability 
timescale. We expect the extrapolation to work better 
in the axial case, where bound-state frequencies have a 
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FIG. 6. (color online) Imaginary part of the axial (left panel) and of the polar S = —1 (right panel) vector modes as a function 
of the BH rotation rate d = J/M^ for I = m — 1 and several values of /i. Although beyond the regime of validity of the first 
order approximation, the modes cross the axis and become unstable precisely when the superradiance condition (|94|l is met. 
As shown in Fig. [T] the first-order results are also in qualitatively good agreement with those obtained at second order. 
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monotonic behavior as functions of a (cf. the left panel 
of Fig. ini). We estimate that close to the onset of the 
superradiant instability the imaginary part of the modes 
should scale as follows: 

Mw, ~ 7s, (am - 2r+/.) {M p)"'^''^''^ , (95) 

where 'yse is a coefficient that depends on S and i. For 
axial modes with i = 1 our numerical data yield 701 ~ 
0.09±0.03, in good agreement with the analytical formula 
in the limit Mfi ^ 1, 

Mojf^' , axial) ^ _L (5m - 2r+fi) {Mfif , (96) 

which is derived in Appcndix[C] From the equation above 
we get 7oi = 1/12 ~ 0.0833. Expression (|96p is compared 
to the numerical results in Fig. [SJ 



Here and in the following, numerical errors are esti- 
mated by comparing the results at first and second order, 
and by taking the maximum deviation between the fit 
and the data. Wc generate numerical data with increas- 
ing accuracy until convergence is reached. In particular, 
in the direct-integration method it is important to retain 
a sufficiently large number of terms in the series expan- 
sion of the boundary condition at infinity, Eq. (j77p . in 
order to achieve sufficient accuracy. This is expected be- 
cause bound-state modes decay exponentially, and there- 
fore it is challenging to perform a numerical integration 
at large distances. 

Wc have applied the same method to compute unstable 
massive scalar modes [60j . whose imaginary part scales 
as in the case of axial vector modes with S = 0. In this 
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case we find ^g^'^''^'' 0.03 ± 0.01, again in reasonable 
agreement with analytical expectations (by setting s = 
in Eq. ^ we get j^l''^''' = 1/48 - 0.0208). Comput- 
ing the modes in the limit Alfi ^ 1 is numerically chal- 
lenging because of the steep scaling of w/ with M ji [cf. 
Eq. (|95p ] and this affects the precision of the fit. Despite 
these numerical difficulties, we are able to recover the cor- 
rect order of magnitude for the instability timescale when 
M/x <C 1. In any case, the observational bounds that we 
will address in the next section depend only mildly on 
the precise value of 75^. For the polar modes we find 
7_ii « 20 ± 10, while extracting the coefficient 711 with 
sufficient precision is numerically challenging even for the 
fundamental mode with ^ = 1, since lji ~ (MfxY^. 

Obviously our method becomes inaccurate when a — >■ 
1, and the fit (|95p can provide at most an order of 
magnitude for the instability timescale in the extremal 
limit. However, in the massive scalar case Eq. (|95p with 
S' = is exact in the M/x <g; 1 limit, for any value 
of a [6^. Furthermore, in the massive vector case, the 
good agreement between the first- and second-order ex- 
pansions suggests that the slow-rotation approximation 
should be quite accurate even for moderately large spins, 
a < 0.7. Therefore it is reasonable to expect that extrap- 
olations of Eq. (PSj) from the slow-rotation limit should at 
least provide the correct order of magnitude (and possi- 
bly a reasonable quantitative estimate) of the instability 
timescale. 

It is tempting to extrapolate our predictions to generic 
values of a and to discuss possible astrophysical impli- 
cations of the superradiant instability for massive vec- 
tor fields around Kerr BHs. Our conclusions rely on 
extrapolation, so they should be confirmed by indepen- 
dent means. Let us assume that the imaginary part of 
the modes is generically described by Eq. (^5]) . where 
we left the coefficient 75^ undetermined, since its precise 
value (whose order of magnitude approximately ranges 
between 0.1 and 10, depending on S and ^) is not crucial 
for the following discussion. From Eq. (pSj) wc find that 
the instability growth timescale has a minimum (i.e., ui 
has a maximum) for 
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(97) 



and the corresponding value of the timescale r 



given by Eq. evaluated at = /x™™. From Eq. 

we also see that (excluding the case am ~ "^r+fj.) this is 

a typical value of the instability timescale. 

If we apply the same procedure to scalar fields {£ = 1, 
S = 0, ''Jqi = 48); we find that the minimum instabil- 
ity corresponds to a = 1, M^™'" — 0.45 and Mu^^^ = 
1.6 X 10~^. Notice that this simple argument relies on an 
(a priori unjustified) extrapolation of analytical expres- 
sions valid for Mji <C 1 to the regime Mfi ~ 1. In the 
scalar case this prediction overestimates the numerical 
results by an order of magnitude: in fact, Refs. [sj, |35| 
found a minimum instability Mujj 1.5 x 10~^ at 
Alfi ^ 0.42 for d ^ 0.99. This disagreement should not 



performs remarkably 
43l |: it over- 



be surprising. Actually, Eq. (|95|) 
well even for values of a close to extremalitv 
estimates the results in Table III in Ref. [35| by only 3% 
when d = 0.7 and by less than 70% when a = 0.99. In 
our view, this agreement is quite impressive. 

Let us apply the same argument to vector fields. From 
Eq. (|95|) we expect the stronger instability when £ = m = 
1. For a = 0.7 this corresponds to 



Mfj. 



™" - 0.187, 
™" - 0.184, 



mill 



0.179. 



5.8711 X 10^1", 
1.7701 X 10-« , 
5.l7_ii X 10-^ 



for S = 1,0,-1, respectively. Our data at second or- 
der are compatible with 701 ~ 0.09 and 7-11 ~ 20. As 
we discussed, we could not extract the coefficient 711 
with sufficient precision. Nevertheless it is not difficult 
to show that the modes with S = 1 are indeed those 
with the smallest (in modulus) imaginary part, i.e. they 
are the least interesting for what concerns the instability. 
This confirms the expectation that polar perturbations 
with S = —1 should have the strongest instability in the 
rapidly rotating limit, as conjectured in Ref. [l^. In the 
Proca case the strongest instability should occur on a 
timescale 



'A'cctor — ^ J 



M{Mfi)- 



7-11 (a - 2^lr+) 



The timescale above must be compared with that for the 
massive scalar case 16(1 for ^ = to = 1, 



^scala 



48M(M/x)- 
d — 2/ir+ 



(99) 



Roughly speaking, the instability timescale against 
vector polar perturbations is of order Tvcctor ~ 
10-27Z,\(M/Me)s. 

In the next section we shall use the results obtained 
from our numerics at second order to discuss some im- 
portant astrophysical consequences of the Proca instabil- 
ity. In this context it is crucial to have a robust estimate 
of the instability timescale in the <^ 1 limit. In 
the axial case our numerical results are supported by the 
analytical formula (pS)) . which provides strong support 
that the fit (^5]) represents the correct behavior when 
Mfi <^ 1. Unfortunately in the polar case we do not 
have the same analytical insight and, as discussed above, 
the small-mass regime is challenging to investigate nu- 
merically. In order to verify the reliability of Eq. (|95p for 
polar perturbations wc have tried several choices of the 
fitting functions in the most relevant case, ^ = to = 1 
and S = —1. It turns out that Eq. (|95p is a conservative 
choice in the polar other fits generically predict 

a stronger instability. Furthermore, by comparing the re- 
sults in Figs. |6] and [7] for the axial and for the polar case, 
it is clear that the polar case exhibits a more complex 
behavior, which is hard to reproduce without some ana- 
lytical insight. In fact our numerical data arc also very 
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well reproduced by the following fitting function: 

MuJi - (a - dsR) ivoiMf^T" + vMMfiT'] , (100) 
where 

-m + y/m^ + 16(A//i)2 ^Mfi 3 
m ^ + ^(^)' 

corresponding to the superradiance threshold [Eq. (|94)) ] 
when ~ /i. The fit (jlOOp has been obtained by im- 
posing a second-order (in a) functional form and ioi ~ 
when 5, = asR, consistently with our data. Finally, the 
functional dependence on fi of the two remaining terms 
has been obtained by fitting the data in the instability 
region and for 0.03 <M^i< 0.1. We found % « -6.5±2, 
771 « 2.1 ± 1, Ko « 6.0 ± 0.1 and ki « 5.0 ± 0.3. While 
the fit ([95l) is physically more appealing [l^, Eq. (|100p 
reproduces our numerical data in the whole instability 
region within a factor of two. It would be interesting to 
better understand the behavior of the polar instability 
in the limit M ^, <^ 1 using different approaches, such 
as a time-evolution analysis or a full numerical evolution 
of the Proca equation. In the following we shall discuss 
how the choice of either Eq. (|95|) or Eq. (|100p affects the 
astrophysical implications of our results. 

VII. ASTROPHYSICAL IMPLICATIONS OF 
THE PROCA INSTABILITY 

Our numerical results in the previous section indicate 
that polar perturbations with S = —1 have the shortest 
instability timescale, as conjectured in [l^. For fixed val- 
ues of a and ^, the instability timescale for polar pertur- 
bations is smaller by two-three orders of magnitude than 
in the axial case. While this conclusion relies on an ex- 
trapolation of calculations that are valid (strictly speak- 
ing) only in the slow-rotation limit, we have shown that 
a similar extrapolation in the scalar case is in remark- 
able quantitative agreement with numerical results that 
do not rely on the slow rotation approximation. Further- 
more, in the Proca case a second-order calculation gives 
solid evidence that our extrapolation is reliable, both 
qualitatively and quantitatively. Therefore it is reason- 
able to expect that extrapolations from the slow-rotation 
limit should at least provide the correct order of magni- 
tude (and possibly a reliable quantitative estimate) of the 
polar instability timescale. Here we explore the tantaliz- 
ing astrophysical implications of the Proca superradiant 
instability. Unless otherwise stated, we shall focus on 
the most relevant modes, those with £ = m = 1, which 
correspond to the stronger instability. 

A. Implications from the existence of supermassive 
black holes 

The most conservative assumption on the final state of 
the instability is that the radiation will extract angular 



momentum from the BH, leaving behind a Kerr metric 
with dimensionless spin parameter below the superradi- 
ant threshold: in other words, the angular frequency of 
the BH horizon must be such that Q,h ^ /i. By this 
argument, even a single reliable supermassive BH spin 
measurement can impose a stringent constraint on the 
allowed mass range of massive vector fields. Bounds on 
the mass of the vector field follow from the requirement 
that astrophysical spinning BHs should be stable, in the 
sense that the timescale (|98p should be larger than some 
observational threshold. 
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FIG. 8. (color online) Bounds on the photon mass m„ = hjj. 
obtained by extrapolating the instability timescale (|98|l for 
S — —1 polar Proca modes of a Kerr BHs up to 5 = 1. 
The region delimited by the curves corresponds to an insta- 
bility timescale r < rnubbic (continuous red line), r < rsaipeter 
(dashed green line) and r < 10* yr (dot-dashed blue line), re- 
spectively. We set 7-11 ~ 20 ± 10 in Eq. (|98p and consider 
a Kerr BH with M = IO^Mq, but the results have a simple 
scaling with 7-11 and with the BH mass (cf. Eqs. H102[) and 
(|103p '). Thin dashed lines bracket our estimated numerical 
errors. 

For isolated BHs we can consider the observational 
threshold to be the age of the Universe, rnubbic = 
1.38 X 10^° yr. A more conservative assumption is to 
include possible spin growth due to mergers with other 
BHs and/or accretion. These effects are in competition 
with superradiant angular momentum extraction. 

The most likely mechanism to produce fast-spinning 
BHf[f| is prolonged accretion [tI]. In this case we can 
compare the superradiance timescale to the shortest 



^ Supermassive BHs with moderate spin (say a ~ 0.7, as in Fairall 
9 [531) could also be produced by comparable-mass BH merg- 
ers [^. The timescale for mergers depends on details of the 
"final parsec problem" , which are poorly known |75|| . The mini- 
mum timescale necessary for BHs to merge is the gravitational- 
radiation timescale, which is of the same order as the Salpeter 
timescale (cf. Eq. (10) in [fl or Eq. (26) in but it is likely 

that the timescale for binary BHs to merge will be dominated by 
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FIG. 9. (color online) Instability timescale (in years) as a function of the vector mass m„ (in eV) for the eight supermassive 
BHs listed in Table 2 of [sH]- The left panel refers to the axial case (5* = 0, ^ = m = 1); the middle panel to the polar case, 
Eq. (|95|) with S — —1, £ — m = 1 and 7_ii ^ 20; the right panel to the polar case, but using the alternative fitting function 
given in Eq. pOOp . The horizontal line corresponds to r = rsaipctor ~ 4.5 x 10^ yr. 



timescale over which accretion could spin up the BH. 
Thin-disk accretion can increase the BH spin from a ~ 
to a « 1 with a corresponding mass increase by a factor 
•\/6 [zl] . If we assume that mass growth occurs via accre- 
tion at the Eddington limit, so that the BH mass grows 
exponentially with e-folding time given by the Salpcter 
timescale rgaipoter = 4.5 x 10^ yr, then the minimum 
timescale for the BH spin to grow from a = to a sa 1 via 
thin-disk accretion is comparable to rsaipotcr- Note that 
accretion at the Eddington limit is a rather conservative 
assumption, since more typically we would expect ac- 
cretion to be sub-Eddington. Furthermore fast-spinning 
BHs are hard to produce in the presence of radiative ef- 
fects [80| or magnctohydrodynamics [8l|, and "chaotic" 
accretion (82j would make large spin parameters even less 
likely. 

For illustration, in Fig. |S] we consider a Kerr BH with 
M = 10^A/q, so that to„ = 10""(M^) eV. We assume 
7_ii sa 20± 10 in Eq. (|M|) (consistently with our numer- 
ical results) and we plot: (i) the superradiant threshold 
(dashed black line), corresponding to the points where 
H « fin and the instability timescale (in the perturba- 
tive treatment) becomes infinite; (ii) the contours cor- 
responding to an instability timescale r = ruubbio (con- 
tinuous red line), r = rsaipotor (dashed green line) and 
T = 10'' yr (dash-dotted blue line), respectively. Each 
curve is shown with the corresponding error bars. 

In the region delimited by the continuous red line the 
instability timescale is shorter than the Hubble time 
(of course, similar considerations apply to the dashed 
green and dash-dotted blue lines). For large a the lower 
branches of the "critical" curves are roughly horizontal. 
By solving ujj^ = Tc in the limit M/i <C 1 this regime 



can be well approximated by the formula 



the dynamical friction timescale, which is typically of the order 
of a few Gyrs for satellites in a Milky Way-type halo [tt} and 
hence larger than the Salpeter timescale. 
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where Tc is the threshold time. If Tc = Tuubbic and in the 
S — —1 polar case, we get 
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Note that, although extracting the values of the fitting 
parameters can be numerically challenging, our results 
depend very mildly on the precise value of 7-11. From 
Eq. (|98p we see that the instability timescale is inversely 
proportional to 7_ii. Therefore the near- horizontal "ob- 
servability threshold" in Fig. [5] that would correspond to 
choosing r = rgaipctor rather than t = : THubbio can be 
obtained multiplying 7_ii in Eq. p03p by the ratio of 
the two timescales, rsaipctcr/THubbic — 3 x 10"'^. Then 
the vector mass bound (jl03|) would increase by a mod- 
est factor of 2.3. What is crucial in obtaining strong 
astrophysical bounds is having a large BH mass, because 
bounds scale like (10^Mq/M)^^^ 

Brcnneman et al. [56j have recently presented a list 
of eight supermassive BH mass and spin estimates. In 
Fig. [9] we show the superradiant instability timescale t 
(in years) as a function of my (in eV) for these supermas- 
sive BHs. The dashed horizontal line denotes Tsaipeter, 
which we (conservatively) assume as the threshold for 
the timescale to be observationally significant. For each 
system we compute instability timescales using the BH 
mass and spin estimates given in Table 2 of [56| (we se- 
lect cither the average spin value quoted in the table, or 
the lower limit on the spin when the observations im- 
ply a > 0.98). The left panel refers to axial modes with 
S = and £ = m = 1, for which we can compute the 
instability timescale analytically: cf. Eq. (|96p and Ap- 
pendix [C] The middle and right panels shows that quali- 
tatively similar results hold for polar modes with S ~ —1 
and £ = m = 1, which exhibit the strongest instability. 
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In the polar case, we computed the instabihty timescale 
using both Eq. ^ (fit I) and Eq. (fit II). The 

results are qualitatively similar and show that, indepen- 
dently of the fitting function, the axial bounds on ruy 
are typically one-two orders of magnitude less stringent 
than the polar bounds. Fig. [9] shows that existing mea- 
surements of supermassive BH spins rule out Proca field 
masses in the whole range lO^'^'^ cV< my < 10~^^ cV. 

Note that arguments based on the superradiant insta- 
bility can only be used to set "exclusion windows" on the 
boson mass, rather then upper bounds [42|, |4^. How- 
ever, observations of BHs in different mass ranges can 
exclude different mass windows. For example, stellar- 
size BHs, M ~ (5 — 20)Mq, can set an exclusion region 
of 10~^^ cV< niy < 10~^ eV. This mass range is relevant 
for so-called "dark photons" coupled to sterile neutri- 
nos [11]. If the identified intermediate-mass BHs were 
confirmed to exist and if the largest known supermas- 
sive BHs with Af ~ 2 X IO^^Mq ^ [HI were confirmed 
to have sufficiently large spin, combined observations of 
spinning BHs can potentially exclude the entire mass 
windows 10-21 eV< to„ < IQ-^ eV [H, |4|. The lat- 
ter is complementary to the bounds set by several other 
experiments and observations that exclude light bosonic 
particle with larger mass [H, HI, [sl] . 

Using current data on supermassive BHs [H^ , the best 
bound comes from Fairall 9 [s^l, for which the polar in- 
stability implies a conservative bound (including mea- 
surement errors) my < IQ-^" eV or niy < 1Q~^^ eV, 
depending on the fitting function (59j . Although the in- 
stability for axial modes is weaker, in this case our re- 
sults are more precise and imply a bound as stringent 
as rUy < 4 X 10-^° eV if we consider r < rsaipotor, or 
my X 10-2° eV if we consider r < rnubbic- Even under our 
very conservative assumptions these results are of great 
significance, since the current best bound on the photon 
mass is nij < IQ-^^ eV [52|,[8^. Combining the exclusion 
window derived above, IQ-^" < niy < 10- eV, with ex- 
isting bounds [52l . l86j , we find the most stringent upper 
bound on the photon mass, < 10-20 eV. 

We remind the reader that our results are also sum- 
marized in Fig. [1] where we show instability windows for 
axial and polar modes together with supermassive BH 
mass and spin estimates from Ref. [56| . 



B. Nonlinear effects, other couplings 

An important ingredient that was not taken into ac- 
count in our study is the nonlinear evolution of the insta- 
bility, that can modify the background geometry. Photon 
self-interactions arc very weak, being suppressed by the 
mass of the electron. Therefore, it is quite likely that 
the outcome of the instability will be a slow and gradual 
drainage of the hole's rotational energy. However, exotic 
massive vector fields may show nonlinearities as those 
studied in Ref. ji^ 113] for the scalar case. 

Another important aspect that should be investigated 



is whether the coupling of accreting matter to massive 
bosons can quench the instability. In principle massive 
photons (unlike hidden U{1) fields, for which the inter- 
action with matter is very small) can couple strongly to 
matter. However it is unlikely that this will significantly 
affect the superradiant instability discussed here, for two 
reasons. The first is that unstable modes are large-scale 
coherent modes whose Compton wavelength is compara- 
ble to (or larger than) the size of the BH's event hori- 
zon, and accretion disks are typically charge-neutral over 
these lengthscales, so any possible coupling with ordinary 
neutral matter is incoherent and most likely inefficient. 
The second is that accretion disks are expected to be 
localized on the equatorial plane, and therefore can af- 
fect at most some (but not all) of the unstable modes. 
The investigation of the superradiant instability in the 
presence of matter requires further work, but the above 
considerations suggest that vacuum estimates should be 
reliable. Spin estimates for slowly accreting BHs (such 
as the BH at the center of our own galaxy) are the most 
reliable, in that they should be less sensitive to details of 
the interaction of vector fields with matter. 



VIII. CONCLUSIONS AND EXTENSIONS 

BH perturbation theory is a powerful tool, with impor- 
tant applications in astrophysics and high-energy physics, 
but its applicability is limited when the perturbation 
equations cannot be separated. We have discussed a gen- 
eral method to circumvent these difficulties in the case of 
slowly rotating BH spacetimes. The method was origi- 
nally developed to study the modes of rotating stars but, 
as we showed, it applies to any slowly rotating BH back- 
ground and to any kind of perturbation. Within this gen- 
eral framework, we discussed two effects induced by ro- 
tation: (i) the Zeeman-like splitting of the eigenfrequen- 
cies, that breaks the degeneracy of modes with different 
azimuthal number m, and (ii) a Laporte-like selection 
rule, in which modes with multipolar index £ are coupled 
to those with multipolar index £±1 and opposite parity 
and to those with £ ±2 and same parity. We have ex- 
tended Kojima's arguments to show that these couplings 
do not affect the QNM frequencies of the spacetime to 
first order in the rotation rate, and we verified this claim 
by numerical integration of the perturbation equations. 

As a first relevant application of this method, we stud- 
ied massive vector (Proca) perturbations on a slowly ro- 
tating Kerr background, i.e. the only interesting class 
of perturbations of Kerr BHs in four dimensions that 
docs not seem to be separable. Using our approach we 
reduced the equations to a set of coupled ODEs, that 
can be solved by a straightforward extension of methods 
valid in the nonrotating case. We have computed the 
corrections to the Proca eigenfrequencies of a Kerr BH 
to first and to second order in the rotation rate. For the 
first time we showed that the Proca bound-state modes 
become unstable when the superradiant condition (|74l) 
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is met. We proved this numerically for the polar modes, 
and both numerically and analytically for the axial modes 
in the limit M/i ^ 1. Unfortunately the method becomes 
intrinsically inaccurate for large values of the rotation 
rate. This prevents us from accurately computing the 
instability timescale close to extremality. By fitting and 
extrapolating our data we confirmed the scaling with to„ 
conjectured in Rcf. [l6j . and we also estimated the order 
of magnitude of the instability timescale. We expect our 
results to be reliable up to a ~ 0.7. Our estimates should 
be confirmed by more accurate methods, but they pro- 
vide strong evidence that measurements of rotating su- 
permassive BH spins set the most stringent upper bounds 
on the photon mass (cf. Fig. [T]). 

Besides their astrophysical implications, our results are 
relevant also for some classes of gravitational theories 
with higher-order curvature terms in the action. For 
example, as proved by Buchdahl [s^, a theory defined 
by £ = (^R + aR[ab]R''"'''^) in the Palatini approach 
is dynamically equivalent to the Einstein-Proca system 
when the vector mass fj,"^ = 3/|a| (see also Ref. jssj). 
The Proca instability of Kerr BHs in this theory can put 
constraints on the parameter a, which is related to the 
nonmetricity of the connection. 
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While limited by the slow-rotation assumption, the 
power of the present method consists in its generality. 
The slow-rotation approximation allows us to study lin- 
ear perturbations of any slowly rotating spacetime even 
when the equations are nonseparable. For instance it 
can be used to study gravitational-electromagnetic per- 
turbations of Kerr-Newman BHs in four dimensions [l2] , 
QNMs of Myers-Perry BHs ^ and BH greybody fac- 
tors in higher-dimensional rotating spacetimes [ill , [s^ . 
Another interesting extension is the stability analysis of 
rotating BHs in asymptotically anti-de Sitter spacetimes. 
This is relevant in the context of the gauge/gravity du- 
ality, as the QNMs of anti-de Sitter BHs in five dimen- 
sions are holographically related to specific correlation 
functions and transport coefficients of the dual bound- 
ary theory 

We note that perturbations of rotating stars have been 
studied up to second order in the Cowling approximation 
(see e.g. Refs. (89l - [9l| ). where metric perturbations are 
neglected and the equations have a structure reminiscent 
of our general Eqs. (|3]) and ([S])- 

Finally we remark that the techniques discussed in this 
paper are not limited to general relativity. They are also 
useful to study gravitational perturbations and the lin- 
ear stability of slowly rotating BH metrics in alterna- 
tive theories of gravity, for example in quadratic grav- 
ity [13, [2^ [2^ . We hope to report on these interesting 
extensions in the near future. 



Appendix A: Coefficients of the perturbation 
equations 

In this Appendix we list the explicit form of the co- 
efficients appearing in the perturbation equations. We 
remark that the coefficients Ai, Di, ,Bf^ ^ of' incor- 
porate terms at zeroth, first and second order in the ro- 
tation parameter a. Therefore, there is no direct corre- 
spondence between these coefficients and the quantities 
at fixed order in a, such as Ai, Ai, Ai, etc., which have 
been introduced in Sec. [Hi In the following, a prime de- 
notes derivative with respect to the radial coordinate r. 

1. Coefficients of the scalar equation 

The explicit form of the coefficients appearing in 
Eq. dHI) is 

Ai = -2 [A'Pr"^ ((4M2 + 2Mr (A - 1 + r^^) + r*^^'^ 
-r^ (A + rV^)) *f + r{r - 2M) {2M% 
-t-r(r - 2M)*^'))] + ^amM^r^uj^ t + h'^M^ [(36A/2 
-H2Afr (A - 11 + r^/?) - (A - 2 + 2m^ + r^/?) 

+r(r - 2M) (lOM^-^ - r{2M + r)*")] , (Al) 
Dt = 2d^M^ [(4A//2 - Ar^ + 2Mr (A - 1 + r^w^)) 

+r(r - 2M) {2M^[ + r{r - 2M)^'[)] . (A2) 



24 



2. Coefficients of tlie Proca perturbation equations 



In this Appendix we list the coefficients appearing in Eqs. (|44p -(|47 p . The coefhcients read 



2r6 



tils (-r^ (A + 2777,2 ^2^2 _ 2) + 32A/2 + 2M7' (A + fi'^r^ ~ lO)) 



2M -r 



(^{2M - r) (^i{e + 1){2M - r) (^-UMu'l^ + r{2M + r)7i"(i) + ir^uju''^^)) ~ 2^"^ r'^ uu^^^ 



£{£ + l){r - 2M) 

-iKr. (SM^ - 2Mr + r') 4,)] + ^ "' \,.(,M - ,') — 

(2M - r) (A + y?r'^) u[^^ + r(r - 2Mfu"\^^ + ((2M - r) (u^g, - + (r - 4M)7i^2) 

r3(2il/- r) 



(A3) 



2Ar4(2il/- 7-)3 



Au[2) ((2M - r) {£^12M + r) + £{2AI + r) + r {2n? + ^?r{2M + r))) + r'^u? (871/^ + 2Mr + r^)) 



+2r772r(r - 2M)2u'(3) + iA?'a;(2A/ - r) (SAf^ - 2Afr + r^) li'^^^ - + l)a;(2A/ - r) (SA/^ _ 6A//r + Sr^) u\^s^ 
2hmAP (^r^uj{r - 2Af)u'(3) + iA (^{2M - r)^^^) + r ((r - 2Af)7i'(i) + 2ira;7i^2)) ) ) 



'"(1) 



{2M - r) {{2M - r)u'(3) + ir'^uju'\^)^ + 71^2) ((2M - r) (A + /i^r^) + r^cj^) - iruj{2M ~ r)7 

r2(r-2A/)2 ' 
a2Af2 (^(2A/ - r) (^-27772rit[3j - Ar(r - 2A.f)77'(2) - A(2Af - r)77[2)) " ^^^^ (^^'^^ - 2A/7' + r^) 77^^^) 



i{£+l){r-2M) 



2 



45777Af2r(rcj77(3) + 7A77[,)) 2r2 (^(2A'/ - r) (r7i'^2) + "(2) - "(3)'( 



A(2M - r) 2Af - r 

4aAf2y^ 4a2777Af^a;77f 



(A5) 



iW = -^!^!:^!^, (A7) 

if=0, (A8) 

r(o) ^ W V IJ_ 12L / Aq^ 

^ £{£ + l)r^2M ~ r) Ar4(2Af-r) ' ^ ' 

nil) _ W w_ (A^()) 

^ £{e + iy{2M-r) Ar2(2A/-r)' ^ ' 

' A(r-2Af) A(r-2A.f) ' ^ ' 
4a2A/3(A77ls + irwTilJ 

— ,(,;v- ' ^^^'^ 

2a2A./2((7--2M)7.-^3)-A4)) 

^ Ar4(2Af-r) ' ^^^^^ 

' = , (A14) 
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(0) 



D 



and 

at 



r6(2Af -r) 
-ira;(10M- r)u[2) 

2M -r 



-{2M - r){6M - Ar)u(i) + r ((2Af - r) ( 



^(3) ^ (2) 



(A15) 



it^2) (-^(^ + l)(2Ai" - r) - 2Afr2cj2) + (2M - r) ((r - 2M)u'(3) - 2afrwM'(i)) + 2iMuj{2M - r)u^^) 



(2A/ - r) (rM'(2) + "(2) - '"(3)^(»')) + 2iA/ra;u^i) 



(A16) 
(A17) 



2kr^{r-2Mf 



-48A'/3u'^3^ + 4Af2r£2u'(2) + 4£Af ^^'(2) + 48AfVu'(3) - AMr^u:'^u\^^ + iAr^uj{2M + r)ufj^-^ 



-4Afr2£^w'(2) - 4£A//r\'(2) - UMr^u'l^) - 3A(r - 2A/)\f2) + 4A/r(r - 2A/r)2u"(3) + r^fu'l^^ + tr 
2amM'^[ruju^,^s—iku^,y?) \ 



u 



(2) 



i{i+l)r^{2M -r) Ar^(2M-r) 



U'l'^u'l^^i - 2^2A//r2u(3) - 2A//r£2u'(2) - 2^Afru'(2) + A(2A.'/ - r)u(2) 

(A18) 



-2A//ra (3) - r(r - 2Af)^w"(3) + M r U(3) - ^ ^^1*^3) + iAr^wu^ij + r^ru (2) + ^r^M (2) 
a2Af2 (£(^ + - 2A/)2 - 2Afr3t^2-) _ 2Af (r - 2Af)2 (ru"l^^ - 3u'^4))) 2dmM 



^(4) 



"(4) 



Ar4(r - 2Af)2 

((2A/ - r) (A + /iV^) + r-^u;2) + (2Af - r) (r(2A/ - r)u"l^-^ - 2A/u'[4)) 



Ar2(2Af-r) 



o = - 



Ar2(2A^-r) ' ^^^^^ 

h^mM^ (^2AA/ra;u[^) - i (^-2A/r2w2u^3j + (r - 2A/) (-8Afu'(3) - r(r - 2M)u"(3) + Aru'(2)) + A(8Af - r)u[2))) 

Ar4(2Af-r) 

6aAf2 (^(2Af - r)u[4) + r ({r ~ 2A/)u'(i) + ira;u|2))) 



2iaM'^u 



\4) 



2Mr^ - r3 



r4(2Af - r) 

m^mM^ (^-2 {6AP - 5A/r + P) u'(4) + 2A/ (A + r^w^) w^^^ + r(r - 2A/)2u"^4) j 

Ar4(2Af -r) 



Ar4(2A-f - r) 



12A'/2u'^3j _ 2MPuj^uI^.^ + 2iKMruju[^^ - 2MrPu'\^) - 2£Mru'\2) + 3A(2A/ - r)u 



(2) 



(4) 



It 



Ar*{2M - r) 



^(4) + 1)(2A/ - r) + 2A//r2a;2) + (2M - r) (^r{2M - r)u"(4) - SAfu' 



(4) 



(A20) 
(A21) 

(A22) 

(A23) 
(A24) 

(A25) 

(A26) 



I 

Appendix B: Proca equation for a slowly rotating first order in 5. Equation ([55]) reads 
Kerr BH 

^2"(4) ;3 "(4) 



6aAf2 



In this Appendix we list the perturbation equations 



l)Qi {Fu[:^^ -irijju\-^^ - Fru 



for a Proca field on a slowly rotating Kerr background to +(-Qi+i [^fujUfl^^ ^"(1) + -^^"^ 



(2) 



(1) 



(1) 



(Bl) 
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where the operator 'D2 is defined in Eq. (|24|) . The Lorenz 
eondition (Eq. ([5^ with 1 = 2) becomes 



2aNPm 
2i aM^u 



F ( ii(2) - "(3) + (2) 



0^(1 



) + T"(3) 



(4) 



(B2) 



Using Eq. (IB2|), Eqs. dM]) with / = 0, 1 and Eq. ^ read 

2aAf2m 



^2"(3) + -2 " 



(2) 



-3iF I U(^) - ru 



0, 



ra;(3u(2) - 2u(3)) 
(B3) 



« 2F f 3M\ 
^2^(2) - 7^ (^1 - — j 



2aM2TO 
6i hhPFuj 

^2"(l)- — (l^"(2) 



"(2) - "(3) 



(^ + l)Q£«(4Y-^Q£+l"(4) 

Fu% 



(3) 



(B4) 



2aM^m 
2aAPF 



e{e + l){2ru}ul^^ -iM[2)) + irFu' 



(3) 



^(£+l)r4 



^ (4) 



(B5) 



The system (|Bip . (jB3|) - (|B5P can be greatly simphfied by 
a systematic use of Eq. (jB2[) to ehminate u^^^^ and U(^)"'^- 
Solving Eq. (|B2p for w^^^ and substituting in Eq. (jBSp - 
(|B4p we get, to first order in a, the equations listed in 
the main text, namely Eqs. (|59|) and (|60|) . 

In order to simplify the axial equation ()B1[) . we con- 
sider Eq. (|B2I) with £ — I and £ + 1 and solve it for u^^J^^ 

and respectively. To first order in a, and making 

also use of Eq. ^ and dSO]), we get Eq. (|6T|) . 



Appendix C: Analytical results for the axial modes 

In this appendix we generalize Detweiler's calcula- 
tion Ig^l of the unstable massive scalar modes of a Kerr 
BH to the case of a massive vector field. We focus on 
the axial sector and, unlike Detweiler, we work in the 
slow-rotation limit up to first order in d. 

As discussed in the main text, the axial vector equa- 
tion (p5|) and the equation for scalar perturbations 
can be written in a unified form, Eq. (j66p . Thus, in 
order to include both scalar and axial vector perturba- 
tions, we shall solve Eq. in the Kerr case, where 
F = B = 1 - 2M/r and vj = 2dAP/r. By defining 



R{r) = ^/r we get 

9 d / 9 dR 
r F— r^F— 
dr \ dr 



r^F £{£ + !) + ^i^r^ 



Admruj 
2Ms2 



i? = 0. (CI) 



For s = 0, the equation above is equivalent to Eq. (7) 
in Ref. [6^ at first order. For s = ±1 and a = it is 
equivalent to Eq. (51) in Ref. [l6| . 

Following Starobinski's method of matching asymp- 
totics [13], we first expand Eq. ((CT|) for r > M: 



d'^jrR) 
dr^- 



2M^2 ^ 



rR = 0. 



(C2) 



The solution of this equation with the correct boundary 
condition at infinity reads 



Roc (x) cx x^e-''/'^U{£ + I - v,2£ + 2,x) , 



(C3) 



where U {p, q, x) is one of the confiuent hypergeomet- 



ric functions ^2. 



and 



v = /koo- The small distance (x <C 1) behavior 
of the solution above reads 

R. 



+ {-ir+'5iy{2£y.i2koor)-'-' . 



(C4) 



where we have defined dv = iy — £— 1 — n and assumed 
,5:/- ^2^+1 < 1 [13. 

On the other hand, Eq. ()C1|) can be solved analytically 
also when r <C niax{£ / lj , £ / fi) . In this limit, by intro- 
ducing a new coordinate z = (r — ?'+ ) /r_|- , the equation 
reduces to 



dz 



Z 



dR 

dz 



[P^ -£{e+l)Z + s^-z]R = 0. (C5) 



We defined Z = z{z + 1) and 

P = -2MkH = -2M{lo - mQ.H) , 



(C6) 



and we neglected 0(5^) terms in P^ . The general so- 
lution of Eq. (|C5[) is a combination of hypergeometric 
functions. By imposing ingoing waves at the horizon, i.e. 
R ~ z"^^ at z ~ 0, we get the general solution with the 
correct boundary condition: 



\2iP , 



Rh{t)o, e-^^^(-l) 

2Fi[-£+-iP + a,l+£ + iP + a,l + 2iP,-z\ , 

where 2Fi{a,b,c, z) is the hypergeometric function [9^ 
and a = — P'^. The large-distance limit of the equa- 
tion above reads 

Ruir) 



(2A/)~'^r[l + 2£] 



r[i + ^ + iP-cr]r[i + ^ + iP + c7] 

(2M)i+^r[-l - 2f] ] 
>riP -(T]T[-£ + iP + a] 



(C7) 
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The key point of the method of matching asymptotics 
is that, when \Mu>\ <^ £ and Mfi <^ £, there exists an 
overlapping region where Eqs. (jC7[) and (jC4[) arc both 



vahd. By equating the coefficients of r and r~ one 
can fix the ratio of the overall factors and the frequency 
of the mode. Finally, we get 



(4fcooAf )^^+ir[-l - 2£]r[2 + n + 2£] ^ 

" ^ r[i + n]r[i + 2^]2r[2 + 2^] ^ 

r[l + iP-cr + £]r[l + iP + (T + £] 
T[iP-(7- £] T[iP + a-£] ' 



(C8) 



prcfactor. Wc find 



(s=±l) 





{i + £y 



(C13) 



which is independent of /i, a and n. The instability 
timescale of the £ = 1 axial vector mode is 4 times shorter 
than that of the scalar mode. Equation (|C12[) is com- 
pared to the numerical results in Fig. [SJ 

On the analytical result in the massive scalar case 



As in the case of massive scalar perturbations, the real 
and imaginary parts of the mode are related to fj. and Siy 
by M 



2 2 2 1 ^^M 

/" -^R = l^ 



£ + n+l 



Siy ( 
i(jji — — 



M \e + n + l 



(C9) 
(CIO) 



In the slow-rotation limit we can further simplify 
Eq. dCl]). To first order in P we get 



iP(4fcooM)2^+ir[2 + 2£ + n] 
r[l + 2£]^T[2 + 2£]^r[l + n] ^ 

s 

r[i + ^ - s]r[i + £ + s]r[i + £]^ J| 



j-l-s-£ 

s-]-£ , 



(Cll) 



When s = we recover the result for the scalar case [60| 
to first order in P (but see the subsection below for a 
discussion). In addition, our results correctly reduce to 
those found in Ref. in the nonrotating limit. We can 
now evaluate the instability timescale. The fundamental 
mode reads 



Mu^r^' = + ^'^■[' + (am - 2r+^) [M^^f . 



Mlo 



48 



(,=2) _ [(2-g)!(2 + .)!] 



(C12) 



885735 



(am — 2r_|-/i) (Mfi) 
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for £ — 1 and £ = 2, respectively. Note that utj ^ /i'^^^^ 
for any value of s, so the only difference arises in the 



Although this is not directly related to the Proca prob- 
lem, we wish to remark that when s = our equa- 
tion ((CTTI) differs from Eq. (28) in Ref. ^ (expanded 
at first order) by a factor 2 for any £ and any n. We be- 
lieve this is either due to an inconsistent limit (explained 
below) or to a typo in Ref. [60| which propagates from 
Eqs. (21)-(24) to Eq. (25), where a factor 1/2 is missing. 
Thus, Eq. (29) in Detweiler should read 7 = fia{fiM)^/48 
for consistency with our Eq. (|C12p . 

In order to understand the source of this discrep- 
ancy we refer the reader to the work by Furuhashi and 
Nambu [o^l, who extended Detweiler's calculation to the 
case of Kerr-Newman BHs. reproducing the results of 
Ref. [gO] in the zero-charge limit. Ref. (sJ] presents a 
detailed derivation of the instability timescale, facilitat- 
ing intermediate comparisons with our own calculation. 
We note that the results in Ref. [U seem to be par- 
tially flawed because the limit defined in Eq. (23) in 
Ref. [13, r(-n)/r(-m) = (-l)"-"m!/n!, is valid only 
when n and m are integers. However, in a rotating back- 
ground the angular eigenvalues are not (in general) in- 
tegers. For example, in the scalar case the separation 
constant \ £{£ + 1) + 0{aMuj) In the Muj < 1 

limit, this can be accounted for by considering £ £ + e, 
with e < 1. When s = 0, the evaluation of Eq. ([CS]) 
above involves limits of the form 



, r(-2^-2e-l) 
lim 

e^o r(-£-e) 



(2£+l)! 



(C14) 



which differs by a factor 2 from Eq. (23) in Ref. [9J] (see 
also the discussion in Ref. [1^). We believe that this is 
the reason for our factor 2 discrepancy with respect to 
Furuhashi and Nambu (and presumably also with respect 
to Detweiler's calculation). Despite numerical difficulties 
in the Af /i <C 1 limit, our numerical calculations for £ = 
m — 1 are in better agreement with a prcfactor of 1/48, 
as computed in Eq. (|C12[) : see also Sec. IVIFl and Fig. 6 
in Ref. M. 
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